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Preface 


The  purpose  of  this  thesis  was  to  develop  a  mechanical  analysis  tool  for  composite 
materials  that  models  the  matrix  and  fibers  separately.  We  intended  this  tool  to  be 
suitable  for  viscoelastic  analysis,  since  the  aerodynamic  heating  anticipated  for  future 
high-speed  aerospace  applications  is  expected  to  cause  either  the  matrix,  the  fibers,  or 
both,  of  the  most  promising  composites  to  behave  viscoelastically.  However,  we  saw 
other  possible  applications  of  this  model,  including  being  able  to  account  for 
delaminations,  incorporation  of  piezoelectric  materials  into  a  composite  panel,  and 
forming  the  basis  of  a  new  type  of  finite  element. 

Once  we  got  started,  we  concluded  that  we  needed  a  rigorous  demonstration  of 
the  validity  of  this  model,  which  required  a  slight  change  in  direction  during  my 
research.  Also,  I  encountered  a  few  road  blocks  and  blind  alleys  along  the  way.  Both  of 
these  kept  me  from  accomplishing  as  much  as  I  had  originally  thought  I  would,  but  I  feel 
that  I  have  established  a  solid  foundation  upon  which  other  research  can  build. 

This  research  was  sponsored  by  Mr.  Jerome  Pearson,  WL/FIBGD,  and  Dr.  T. 
Nicholas,  WL/MLLN. 

I  want  to  thank  my  thesis  advisor,  LtCol.  R.L.  Bagley,  for  helping  me  find  the 
way  around  those  roadblocks  and  blind  alleys,  and  for  putting-up  with  some  of  my  dumb 
mistakes.  I  also  want  to  thank  Drs.  P.J.  Torvik  and  A.N.  Palazoito,  as  members  of  my 
thesis  committee,  for  their  contributions. 


Alan  L.  Lesmerises 
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Abstract 


^  A  new  model  for  composites  represents  a  composite  panel  as  an  alternating  stack  of 
isotropic  matrix  and  orthotropic  fiber-dominated  layers.  An  energy  approach  results  in 
a  set  of  partial  differential  equations  and  boundary  conditions  where  the  mechanical 
properties  of  the  matrix  and  fibers  appear  separately.  Two  solution  methods  are 
developed.  The  first,  called  Strata  Theory"  is  developed  for  general  applications,  while 

A  ^ 

a  second,  called  the  "Exact  Stratified  Plate",  is  used  as  a  benchmark  for  comparison  to 
Strata  Theory.  Both  solutions  assume  zero  gradients  of  displacement  through  the 
thickness  of  a  fiber  layer.  They  are  compared  for  an  isotropic  material  and  a  two-ply 
composite  with  cross-plies,  using  simply-supported  cases  of  a  square  plate  ano  cylindrical 
bending. 
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DEVELOPMENT  OF  A  STRATIFIED  PLATE  MODEL 


FOR  COMPOSITE  PANELS 


1.  Introduction 


Background 

A  composite  material  typically  consists  of  a  relatively  soft  matrix  material  that 
contains  stiffer  fibers.  The  length  of  the  fibers  used  may  vary,  depending  on  the 
application,  from  very  short  fibers  (i.e.  particles)  to  long  continuous  strands.  These  fibers 
endow  the  composite  material  with  increased  mechanical  stiffness.  The  increase  in 
stiffness  depends  primarily  on  the  fiber  density,  properties,  and  orientation. 

When  the  composite  is  formed  into  a  panel,  it  is  usually  built-up  as  a  series  of 
layers  of  parallel  fibers  known  as  plies.  If  each  fiber  J.T  y  is  nodeled  as  an  orthotropic 
plate,  the  panel  as  a  whole  could  be  modeled  as  a  set  of  tiiese  orthotropic  plate  stacked 
upon  one  another.  In  fact,  a  commonly  used  method  for  the  analysis  of  composite  panels, 
known  as  Classic  Laminated  Plate  Theory  (CPT),  models  laminated  composites  this  way. 
In  CPT,  the  elastic  properties  of  the  matrix  are  "averaged"  with  those  of  the  fibers,  and 
the  ability  to  separately  model  the  matrix  and  the  fiber  materials  is  lost. 

At  high  temperatures,  many  of  the  materials  used  for  the  matrix  and  fibers  in 
today’s  composites  can  become  viscoelastic.  In  such  a  case,  there  is  no  reason  to  expect 
that  the  fibers  and  matrix  will  become  viscoelastic  at  the  same  time,  or  that  the  fibers  and 
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matrix  would  maintain  the  same  elastic  relationships  once  they  do  become  viscoelastic. 
Therefore,  it  is  essential  that  any  type  of  vibration  analysis  where  viscoelasticity  will  be 
considered  is  able  to  deal  with  the  matrix  and  fiber  elasticity  terms  separately. 

Scope 

There  is  a  need  for  an  alternate  method  of  modeling  composites,  such  that  the 
matrix  and  fiber  elasticity  terms  appear  separately  in  the  governing  equations. 
Furthermore,  a  formulation  suitable  for  viscoelastic  stress  and  vibration  analyses  is 
needed. 

This  thesis  details  the  development  of  a  new  model  for  the  analysis  of  laminated 
composites  called  the  Stratified  Plate  Model.  A  composite  is  modeled  as  an  alternating 
stack  of  isotropic  matrix  strata  and  orthotropic  fiber  strata.  As  a  result,  the  matrix  and 
fiber  properties  have  an  identifiable  and  distinct  contribution  to  the  overall  mechanical 
characteristics  of  the  composite. 

A  new  analysis  technique,  called  Strata  Theory,  arises  from  the  Stratified  Plate 
Model  and  is  intended  for  general  applications.  Strata  Theory  will  be  validated  by 
comparing  it  with  a  more  accurate  formulation  called  the  "Exact  Stratified  Plate". 

Idealization  of  Fiber  Layers 

The  Stratified  Plate  Model  treats  each  fiber  ply  somewhat  differently  than  the 
CPT  model.  First,  each  fiber  is  idealized  [Figure  1.]  with  a  square  cross-section  of  equal 
cross-sectional  area.  The  thickness  of  a  fiber  stratum  hf  is  then  given  as 
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where  Df  is  the  average  diameter  of  the  actual  fibers.  The  shape  of  the  fibers  now  define 
three  distinct  regions  within  the  ply:  a  fiber-dominated  stratum  sandwiched  by  two  thin 
matrix  strata  above  and  below  [Figure  2.].  The  matrix  strata  consist  of  purely 
homogeneous  matrix  material,  while  the  fiber  strata  contains  the  idealized  fibers  and  the 
matrix  material  between  the  fibers  within  the  ply. 


Composite  Panel 


Stratified  Plate  CPT  Model 


Figure  1.  Comparison  between  Stratified  Plate  and  CPT  models  of  a  composite  panel 

Notice  that  when  the  fibers  are  modeled  this  way,  the  idealized  fibers  span  the 
entire  thickness  of  the  fiber  stratum.  Since  the  fiber’s  elastic  modulus  is  typically  at  least 
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Figure  2.  Strata  Theory  idealization  of  fibers 


an  order  of  magnitude  larger  than  that  of  the  matrix,  it  seems  reasonable  to  assume  that 
the  z-direction  displacement  gradients  through  the  fiber  strata  will  be  much  smaller  than 
those  through  the  matrix.  Indeed,  this  assumption  is  a  key  feature  of  the  Stratified  Plate 
Model.  This  will  become  important  in  the  development  of  the  governing  partial 
differential  equations  (PDEs),  as  it  will  allow  all  the  z-direction  displacement  derivatives 
through  the  fiber  strata  to  be  set  to  zero. 

Development  and  Validation 

The  governing  differential  equations  and  associated  boundary  conditions  are 
derived  using  the  principle  of  minimum  potential  energy.  The  Strata  Theory  model  is 
developed  using  an  assumed  displacement  field.  The  Exact  Stratified  Plate  model,  on  the 
other  hand,  finds  the  z-dependent  portion  of  the  displacement  field  through  an  eigenvalue 
formulation  once  a  set  of  x-  and  y-dependent  functions  are  found  that  satisfy  the 
boundary  conditions. 

Exact  orthotropic  plate  solutions  have  been  found  for  composites  in  cylindrical 
bending  and  for  rectangular  bidirectional  sandwich  composite  plates  [8:398-411;  9:20- 
34].  These  solutions  were  used  for  comparison  with  CPT.  Both  CPT  and  these  exact 
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orthotropic  plate  solutions  model  a  composite  as  a  set  of  stacked  orthotropic  plates.  This 
is  considerably  different  from  the  way  the  Stratified  Plate  Model  represents  a  composite. 
These  two  different  representations  of  a  composite  may  predict  radically  different 
behaviors  for  the  composites  they  model.  As  a  result,  it  was  felt  that  using  results 
obtained  from  either  CPT  or  the  exact  orthotropic  plate  methods  as  a  basis  of  comparison 
for  Strata  Theory  would  not  be  appropriate.  Strata  Theory  needed  a  different  standard 
against  which  it  could  be  compared,  hence  the  development  of  the  Exact  Stratified  Plate 
solution. 


Figure  3.  Relationships  between  the  various  models  for  composites 
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II.  The  "Strata  Theory 


Basic  Concepts  and  Assumptions 

The  Strata  Theory,  hereafter  referred  to  as  ST,  is  based  on  the  Stratified  Plate 
Model  described  above.  The  Stratified  Plate  Model  assumes  that  the  z-direction 
displacement  gradients  through  the  fibers  are  zero.  This  means  that  et  will  be  zero 
through  the  fibers.  The  ST  model  makes  the  additional  assumption  that  e2  is  zero  through 
the  matrix  strata  as  well.  Since  ST  was  intended  to  form  the  basis  of  a  new  vibration 
analysis  tool  and  not  for  detailed  stress  and  strain  estimation,  it  seems  to  be  a  reasonable 
approximation  that  the  contributions  of  e2  from  the  matrix  strata  will  be  negligible. 

Another  approximation  is  that  ST  assumes  that  the  displacement  field  can  be 
represented  as  a  polynomial  in  z,  leaving  the  x  and  y  dependencies  as  unknowns  that 
need  to  be  found. 

The  presence  of  fibers  in  a  composite  endow  the  composite  with  complex 
mechanical  properties  such  that  coupling  may  occur  between  the  fiber  plies.  When  there 
is  coupling,  it  means  that  plane  sections  neither  remain  plane  nor  perpendicular  to  the 
midplane.  In  other  words,  what  would  start  as  a  vertical  line  through  the  thickness  of  an 
unloaded,  undeflected  composite  panel  will  become  a  complex  3-dimensional  curve  when 
the  panel  is  deformed.  To  represent  this  complex  3-D  curve,  the  following  expressions 
for  the  displacements  u,  v,  and  w  will  be  used  for  displacements  through  the  fibers 
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(2) 


N 

fay) 

i- 0 
N 

vf(x,yji)  =  ^zfCfay) 
i- o 

*yCw)  =  H<x,y) 

while  the  displacement  field  through  the  matrix  will  use  the  following  expressions  for  u, 
v,  and  w. 

N 

um(x,y*)  =  £(«m+Pmzy  *,(*,y) 

i- 0 

N  (3) 

vm(w)  =  Ci(*»y) 

i-0 

wM(x^)  =  w(x,y) 

The  term  zf  that  appears  in  the  fiber  displacement  terms  above  is  the  z  coordinate  of  the 
midplane  of  a  fiber  stratum  referenced  to  the  midplane  of  the  entire  stack  [Figure  4.]. 
The  values  of  z  that  are  used  in  the  matrix  displacement  field  expressions  must  be  a  z 
coordinate  that  lies  within  a  matrix  stratum.  Also,  the  w,  \J/it  and  £  terms  that  appear  in 
the  assumed  displacement  field  expressions  are  the  same  for  both  the  matrix  and  fibers. 

The  complex  3-D  curve  described  above  can  be  decomposed  by  taking  projections 
onto  the  x-z  and  y-z  planes.  The  in-plane  displacements  u  and  v  are  polynomials  in  z  that 
approximate  these  projections,  and  the  and  £  terms  are  the  x  and  y  dependent 
coefficients  of  these  polynomials  in  z.  It  should  be  noted  that  the  number  "N",  the  order 
of  the  assumed  polynomial,  is  not  directly  related  to  the  number  of  fiber  plies.  However, 
the  more  twisting  and  warping  anticipated,  the  larger  N  should  be  to  better  approximate 
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this  deformation.  Ultimately,  N  is  dictated  by  the  accuracy  required,  and  the  complexity 


of  the  problem. 

4Z 


Figure  4.  Definition  of  the  fiber  strata  parameter  z,  with  respect  to  the  z  axis  origin 

The  a’s  and  /S’s  that  appear  in  the  matrix  displacement  functions  are  mapping 
terms.  Since  the  z  gradients  through  the  fibers  are  assumed  to  be  zero,  the  x  and  y 
displacements  are  constant  through  a  given  fiber  stratum.  This  implies  that  there  must 
be  continuity  of  in-plane  displacements  (u  and  v)  from  the  bottom  of  one  matrix  stratum 
to  the  top  of  the  next  lower  matrix  stratum.  This  continuity  of  displacement  is  established 
by  a  linear  mapping  of  the  z  coordinate  of  the  lower  surface  of  the  ith  matrix  stratum  and 
the  z  coordinate  of  the  upper  surface  of  the  (i+l)th  matrix  stratum  to  the  z  coordinate 
of  the  midplane  of  the  intervening  fiber  stratum,  denoted  by  Zf. 


a 


°v  P-S  =  *f 


(4) 


Each  matrix  stratum  has  its  own  a  and  0.  Finding  a  and  0  for  a  given  matrix 
stratum  is  a  simple  linear  problem  of  2  equations  and  2  unknowns.  The  first  equation 
maps  the  z  coordinate  of  the  upper  surface  of  the  matrix  stratum  to  the  midplanc  of  the 
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fiber  stratum  immediately  above;  the  second  equation  maps  the  z  coordinate  of  the  lower 
surface  of  the  matrix  stratum  to  the  midplane  of  the  fiber  stratum  immediately  below. 
For  a  matrix  stratum  at  the  top  or  bottom  of  the  composite,  the  free  surface  coordinate 
(h/2  or  -h/2)  is  used  in  place  of  zf. 


Vmfa  =  Zfu  (uPP*r  fiber  tyer) 
arn,+  Vmfa  =  Zfi  Oower  fiber  layer) 


(5) 


When  all  the  fiber  plies  of  the  original  composite  being  modeled  are  of  equal  thickness, 
0  is  the  same  for  all  the  matrix  strata 


where  tp  is  the  thickness  of  the  original  fiber  ply.  Note  that  this  is  not  a  requirement  for 
the  ST  approach.  However,  for  the  sake  of  simplicity,  this  assumption  will  be  used  for 
the  remainder  of  the  ST  development. 

These  assumed  displacement  field  expressions  can  now  be  substituted  into  the 
expressions  for  strain,  which  will  later  be  used  to  determine  the  strain  energy 

0 
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Stress  and  strain  for  the  matrix  strata  are  related  through  the  expression 
{a}=[K]{e},  where  the  elasticity  matrix  [K]  for  the  matrix  material  is  given  by 
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or  in  simplified  notation 


i-'ml  Vml 

Vml 

0 

0 

0 

< 

8. 

— 

i 

< 

8. 

Vml 

0 

0 

0 

Vml  V*1  1~Vml 

0 

0 

0 

0  0 

0 

l-2vw 

2 

0 

0 

0  0 

0 

0 

1-2V* 

2 

0 

0  0 

o 

0 

0 

l-2v 

V 

2 

^  ** 

0 

0 

0 

K  K  K 

0 

0 

0 

K  ^ 

0 

0 

0 

0  0  0 

*3 

0 

0 

0  0  0 

0  j 

0 

0  0  0 

0 

0  J 

K, 

*ml  J 

(9) 


(10) 


10 


where 


Kt  = 


(l-2vj(l+vj 


Y 

(1-2vJK1*vJ 


JY  _  ini 

3-_  2(T^J 


(11) 


As  discussed  previously,  the  fiber  strata  are  modeled  as  orthotropic  layers.  The 
elastic  constants  for  such  a  layer  are  derived  using  a  "rule  of  mixtures".  For  this 
derivation,  the  following  fiber  elasticity  terms  will  be  defined 


Sx  =  vol, 


Efd-Vf) 

- L - L -  +  \ol 

f  (1-2V/X1+V/) 


g,a-y,) 

(1-2V.X1+V.) 


s  =  i  *.g-o 

2  vo/»  (l-2v.Xl*v.) 


S3  =  vol, - 22 -  +  vol  - ^ - 

3  /(l-2v/)(l+v/)  m(l-2vm)(Uvm) 

S  =  £»v» 

4  "  (1-2V.X1+V.) 

I  volm 

S6  =  volj  Gf  +  w/,,  Gm 

*7  =  <?*  / 


(12) 


The  terms  volf  and  vo^  are  the  volume  fractions  of  fibers  and  matrix  in  a  given  fiber 
stratum.  This  must  not  be  confused  with  the  global  volume  fractions  of  fibers  and  matrix, 
but  there  is  a  direct  relationship  between  the  two. 

The  elastic  stiffness  for  each  fiber  stratum  in  its  principal  axis  directions  will  have 
the  form  [2:35] 
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where  the  "1"  direction  is  parallel  to  the  fibers,  and  the  "2"  direction  is  perpendicular 
to  the  fibers  in  the  plane  of  the  fiber  stratum.  The  coordinate  transformation  matrix  [0] 
for  an  orthotropic  material  is  given  by  [2:48-49] 
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This  transformation  can  be  used  for  stresses  as  well.  This  is  now  applied  to  the  principal 
axis  elasticity  matrix  by  means  of  the  similarity  transformation 

(«),  *  [e)[s][9]‘H  <15) 


where  [0]  is  the  rotation  matrix  above.  This  produces  the  general  constitutive  relation 
[2:34] 
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Figure  5.  Definition  of  fiber  orientation  angle  6f 


Development  of  the  PDEs 

The  development  is  based  on  the  principle  of  minimum  potential  energy,  given 


by 


// 


♦  ♦ 

H  *"1  F  *** 

£  f  dz-j^  f 6IU  dz  +p(x,y)  t>w 
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drdy  =  0 


(17) 


The  expressions  for  the  strain  energy  for  the  matrix  and  fiber  strata  are  given  by 
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(18) 


n„,  = 


where  F  =  the  number  of  fiber  strata,  and  M  =  the  number  of  matrix  strata.  The  strain 
energy  is  derived  using  the  assumed  displacement  field  expressions,  along  with  the  [K] 
and  [S]  matrices  given  above.  The  work  and  strain  energy  expressions  are  combined,  and 
the  principle  of  minimum  potential  energy  is  applied. 

Since  the  z  dependence  has  already  been  assumed,  the  z  derivatives  that  arise 
from  the  expressions  for  strain  can  be  evaluated  directly,  and  integration  with  respect  to 
z  can  be  performed.  However,  integration  by  parts  is  necessary  to  produce  the  boundary 
conditions,  and  to  eliminate  the  x  and  y  derivatives  of f*-  *  variations  on  w,  and 
The  governing  PDEs  look  like  higher  order  piate  theory  with  coefficients  that  depend  on 
material  (matrix  or  fiber)  properties.  The  PDEs  that  result  are 
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and 
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i*0 


where 


u  z-' 


There  are  2N+2  sets  of  boundary  conditions  that  result  from  this  energy 
derivation,  and  they  are 
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111.  The  " Exact  Stratified  Plate 


Basic  Concepts  and  Assumptions 

The  Exact  Stratified  Plate  model,  hereafter  referred  to  as  ESP,  is  intended  to  be 
a  standard  against  which  the  results  of  ST  can  be  compared.  Both  ST  and  ESP  are  based 
on  the  same  Stratified  Plate  Model,  where  it  is  assumed  that  z-direction  displacement 
gradients  are  zero  through  the  fiber  strata.  Further,  ESP  is  somewhat  similar  to  ST  in 
that  the  differential  equations  are  developed  using  the  principle  of  minimum  potential 
energy.  However,  that  is  where  the  similarity  ends. 

There  are  two  primary  differences  between  ESP  and  ST  methods.  The  first  is  that 
ESP  allows  ez  to  be  non-zero  through  the  matrix.  Recall  that  one  of  the  first  assumptions 
of  ST  was  that  ez  could  be  neglected  in  both  the  fiber  and  matrix  strata.  This 
approximation  seems  reasonable,  but  lacks  analytical  justification.  By  allowing  ez  to  be 
nonzero,  ESP  will  be  useful  for  estimating  the  effects  of  assuming  ez=0. 

The  second  major  difference  between  ESP  and  ST  is  that  ESP  is  restricted  to 
modeling  only  simply-supported  composites  with  cross-plies.  This  limitation  arises 
because  of  the  choice  of  the  assumed  displacement  field.  For  this  analysis,  ESP  assumes 
that  the  form  of  the  x  and  y  dependence  of  the  displacement  field  is  the  same  as  the  ST 
displacement  field  expressions  used  for  the  simply-supported  example  problems.  This  is 
not  really  a  drawback,  it  simply  means  that  ESP  and  ST  can  only  be  directly  compared 
for  a  simply-supported  plate. 

Another  consequence  of  using  this  form  of  the  assumed  x  and  y  dependence  is 
that  it  cannot  accommodate  angle  plies.  In  developing  the  governing  differential 
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equations,  all  coupling  terms  that  arise  from  having  angle  plies  would  not  have  the  same 
x  and  y  dependence  as  the  rest  of  the  terms.  If  these  coupling  terms  were  nonzero,  the 
x  and  y  dependent  functions  would  not  result  in  the  set  of  ODEs  that  result  from  this 
analysis.  Therefore,  only  a  composite  with  cross-plies  can  be  modeled  using  the  ESP 
approach. 

For  ST,  it  was  assumed  that  the  z  dependence  of  u  and  v  could  be  described  by 
a  polynomial  in  z.  However,  the  ESP  solution  leaves  the  form  of  the  z  dependence 
through  the  matrix  strata  as  unknown  functions  that  are  determined  by  the  system  of 
ODE’s  that  result.  Since  this  approach  lets  the  problem  dictate  what  the  z  dependence 
should  be,  ESP  can  be  used  to  gauge  how  well  the  polynomial  used  in  ST  represents  the 
more  complex  stress,  strain,  and  displacement  behavior  that  a  composite  would 
experience. 


Development  of  the  ODEs 

All  strain  terms  are  treated  as  being  nonzero  for  the  matrix,  but  strains  through 
the  fibers  are  assumed  to  have  no  z  dependence.  Therefore,  the  expressions  for  the 
strains  become 
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The  ESP  solution  is  developed  using  the  principle  of  minimum  potential  energy, 


18 


just  as  with  ST.  The  strain  energy  is  initially  developed  using  the  expressions  for  the 
strains  shown  above.  Integration-by-parts  is  performed  to  eliminate  the  derivatives  of  the 
variational  terms.  This  results  in  a  set  of  PDEs  and  boundary  conditions.  The  integration- 
by-parts  produces  a  set  of  stress/displacement  boundary  condition  expressions  for  the  4 
"edges"  of  the  composite  plate,  as  well  as  for  the  horizontal  (x-y  plane)  surfaces  of  the 
composite.  The  resulting  expression  for  the  energy  is 


“  £  6Umi  ~  £  bufk  +  ffp(x,y)6wdxdy  =  0 


where  p(x,y)  is  applied  at  the  top  and/or  bottom  surface  (±h/2), 

6  aO  J<vsy+wJ]6u 

6U-‘  =  ~fff  dzdydx 

bzml  azml  b 

+  If  [*i  +  / /  *>J*,+*j*»  ** 

o  _ -  o - 

zml  *ml 

♦  ♦ 

azml  J  b^ml 

+  If  [^+V,+^]6v|0^  +  I!  *3 

o.-  o  - 

«  b  t  bz ml 

+  If  *  If  K3Ju^)bw  dz4y 

0  0  0j^ 

a  b  a  b 

I  ♦  I  ♦ 

+  JJ  KsJtUj+wJbuV^dydx  +  JJ  ^Jv^+wJbvV^dydx 


If  KiJyswJ*w |0** 


19 


b  a 

du 

'T* 

ii 

•o 

+ 

V-»+V»+(VV,» 

5v 

00 

+ 

'V,.»+V’®]®w 

dydx 


/VV^V->]S“^  *  (32) 

0  0 

a  b 

NVA^K*  +  /VV“^v->)Svlo^ 

0  0 

*  a 

/W^8w|!4’  *  /W/*!** 


The  term  hft  represents  the  thickness  of  the  k*  fiber  stratum,  which  arises  from  the 
integration  with  respect  to  z  through  a  fiber  stratum. 

This  energy  expression  results  in  three  "field  equations",  a  set  of  coupled  PDEs, 
for  each  matrix  stratum.  The  remaining  terms  form  the  boundary  conditions  for  this  set 
of  PDEs  (see  discussion  below). 

As  described  above,  the  x  and  y  portions  of  the  displacements  u,  v,  and  w 
through  the  matrix  are  assumed  to  have  the  same  form  as  the  ST  displacement  field 
expressions  for  a  simply-supported  composite  with  cross-plies.  The  assumed  displacement 
field  and  loading  expressions  are  given  by 

vCw) 
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w(x,y,z) 

P(x,y) 


U(z)cos(px)sm(qy) 
P(z)  sin  (px)  cos(qy) 
JP(z)sin(px)sin(qy) 
Psin(px)sin(0») 


where  p  =  r/a  and  q=x/b.  Since  the  displacements  u,  v,  and  w  for  the  fibers  are  assumed 
to  have  no  z  dependence,  U(z),  V(z),  and  W(z)  are  constant  through  the  fiber  strata. 


20 


The  form  of  the  surface  load  "P"  dictates  this  choice  of  the  assumed  displacement 
field.  These  displacement  field  expressions  and  the  subsequent  differential  equations  can 
be  generalized  for  more  complex  surface  loads  using  Fourier  analysis. 

When  these  displacement  field  expressions  are  substituted  into  the  PDEs,  the  x 
and  y  derivatives  can  be  evaluated  directly,  and  a  set  of  ODEs  for  U(z),  V(z),  and  W(z) 
result.  The  functions  that  are  used  for  the  variational  analysis  are  the  U(z),  V(z),  and 
W(z)  of  the  assumed  displacement  field,  and  the  x-  and  y-dependent  functions  can  be 
factored  out  to  produce  the  following  3  second-order  coupled  ODE’s  for  each  matrix 
stratum  T. 

-  o 

twiwwv  - 0 


The  solution  of  this  set  of  ODEs  is  obtained  by  an  associated  Eigenvalue  formulation. 

The  field  equations  for  each  matrix  stratum  "1"  can  be  manipulated  into  matrix 
form  to  become 
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where  K4=K2+K3.  These  can  be  recast  as  a  set  of  six  first-order  coupled  ODEs  to  get 


0  0  0  A^  0 

0  0  0  0  A , 

0  0  0  0  0 

A^  0  0  0  0 

0  A,  0  0  0 

0  0  Klmi  -pK4mi  -qK4mi 

A,  0  0  0 

0  A,  0  0 

0  0  A.  0 

lml 

0  0  0  p'K^+q1*^ 

0  0  0  pqK4mi 

0  0  0  0 


If  the  solutions  for  U,  V,  and  W  are  assumed  to  have  the  form  of  U^,  the  governing 
equations  can  be  expressed  as  a  system  of  6  coupled  algebraic  equations 


22 


Rewriting  this  in  state-space  form,  this  becomes 


where 


and 


[[<?]  +  M*]]{0O}  -  {0} 


(38) 


(39) 


The  eigenvalues  of  this  system  are  two  sets  of  three  repeated  roots.  The  two  values  of 
the  roots  are 


Note  that  these  roots  are  independent  of  the  material  stiffness  coefficients. 

The  repeated  roots  means  that  an  ordinary  eigenvector  solution  technique  cannot 
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be  employed  here.  However,  Chen  [1:53-57]  outlines  a  more  general  method  for 
generating  an  eigenvector  matrix. 

The  basis  of  this  method  is  as  follows:  let  [A]  be  an  nxn  matrix  with  n 
eigenvalues  (they  do  not  have  to  be  distinct),  and  f(X)  be  a  function  defined  on  the 
spectrum  of  [A].  If  a  polynomial  g(X)  is  constructed  such  that  g(X)  and  f(X)  have  the 
same  values  on  the  spectrum  of  [A],  then  f([A])  can  be  calculated  using  the  polynomial 
defined  by  g([A]).  When  the  n  eigenvalues  are  distinct,  f(X)  and  g(X)  must  be  equated 
to  solve  for  the  coefficients  in  the  polynomial  g(X).  When  there  are  repeated  roots, 
derivatives  of  f(X)  and  g(X)  must  be  equated  as  well. 

For  the  eigenvalue  problem  given  by  ESP,  this  approach  starts  by  defining  the 
following  terms 

g(z,i)  -  «„'-)*«  ,te)  X  ♦  «j(j)  X2  *  a, ft) X3  +  a4(z)  I*  ♦  «,(z)  X! 
f(z  X)  -  ** 

where  the  af’s  are  unknown  functions  of  z,  and  X  is  an  eigenvalue  as  defined  above. 

Since  there  are  repeated  roots,  the  derivatives  of  f(z,X)  and  g(z,X)  are  equated 
such  that 

HzX)  *  *<U)  ,  4-  /fc*>  *  4  •  o'*  -  —,g&»  <«> 

dk  dk  dk1  dX2 

Next,  these  expressions  are  manipulated  to  solve  for  the  a^s.  For  this  problem,  when 
there  are  two  sets  of  triple  roots,  the  a/s  are  found  by  the  relation 
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Finally,  exponentiation  of  a  matrix  [Q]  can  be  obtained  by  the  relation 

« lQ]l  =  « 0[I]  +«,[(?]  +  a2[Q ]2  +  a3[<?]3  +  a4[Q ]4  +  a3[<?]5  (45) 


Assembling  the  Boundary  Conditions 

The  PDEs  are  solved  subject  to  boundary  conditions  along  the  four  edges  (at 
x=0,  x=a,  y=0,  and  y=b),  at  the  top  and  bottom  horizontal  faces,  as  well  as  at  each 
fiber  stratum.  At  the  fiber  strata,  these  internal  boundary  conditions  represent  stresses 
r„,  ryI,  and  ot  at  the  top  and  bottom  surfaces  of  each  matrix  stratum.  External  loads  (P) 
applied  at  the  upper  and  lower  surfaces  of  the  composite  are  also  incorporated  into  these 
boundary  conditions. 

When  the  x  and  y  portions  of  the  displacements  u,  v,  and  w  through  the  matrix 
are  assumed  to  have  the  same  form  as  the  ST  displacement  field  expressions  for  a 
simply-supported  composite  with  cross-plies,  all  the  edge  boundary  conditions  are 
satisfied  exactly. 

The  assumption  of  no  z-gradients  of  displacement  through  the  fiber  strata  also 
implies  a  set  of  continuity  boundary  conditions.  The  values  of  u,  v,  or  w  are  constant 
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through  a  fiber  stratum,  therefore  values  of  U,  V,  and  W  at  the  bottom  of  the  matrix 
stratum  just  above  this  fiber  stratum  must  be  equal  to  the  corresponding  values  at  the  top 
of  the  matrix  stratum  just  below  the  fiber  stratum. 

Since  U,  V,  and  W  are  constant  through  each  fiber  stratum,  integration  through 
the  fibers  can  be  performed  directly.  The  results  of  the  integration  through  a  fiber 
stratum  are  combined  with  the  boundary  condition  terms  of  the  matrix  strata  immediately 
above  and  below  this  fiber  stratum.  In  other  words,  the  sum  of  all  terms  that  are 
evaluated  at  the  i*  fiber  stratum  (the  terms  at  the  bottom  of  the  upper  matrix  stratum,  the 
fiber  stratum  itself,  and  the  top  of  the  matrix  stratum  below)  must  sum  to  zero.  This 
forms  a  set  of  internal  stress  boundary  conditions  between  the  matrix  strata. 

For  M  matrix  strata,  there  will  be  3(M-1)  internal  stress  boundary  conditions,  and 
the  presence  (or  lack  of)  external  loads  at  the  upper  and  lower  surfaces  of  the  composite, 
i.e.  at  z=h/2  or  z=-h/2,  specify  a  total  of  6  more  stress  boundary  conditions. 
Additionally,  it  is  necessary  to  specify  3(M-1)  displacement  (continuity)  boundary 
conditions  to  tie  the  bottom  surface  of  one  matrix  stratum  to  the  top  of  the  next  matrix 
stratum  below  it.  In  all,  there  will  be  a  total  of  6M  boundary  conditions  (both  stress  and 
continuity)  specified  for  M  matrix  strata  that  will  be  used  to  solve  for  the  functions  U, 
V,  and  W. 

The  exponential  expressions  for  U,  V,  and  W  describes  how  U,  V,  and  W  vary 
with  z,  subject  to  the  boundary  conditions  specified  in  Equations  (31)  and  (32).  The 
boundary  conditions  are  assembled  into  a  series  of  stress  balance  equations  for  r„,  ryi, 
and  <t,.  Equilibrium  is  established  at  the  upper  surface,  at  each  fiber  stratum,  and  finally 
at  the  lower  surface.  For  a  composite  with  two  fiber  plies,  this  would  be  given  as 
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The  submatrices  Ml,  and  M2,  are  the  boundary  conditions  given  in  Equation  (31) 
expressed  in  matrix  form.  Submatrix  M2&  incorporates  the  interstrata  boundary 
conditions  imposed  by  the  fiber  strata  given  in  Equation  (32).  Continuity  of  displacement 
across  the  fiber  strata  is  established  by  the  3x3  identity  submatrices,  denoted  by  "I" 
above. 


Determining  the  solution  for  U(z),  V(z),  and  W(z)  depends  on  finding  a  U„,  V0, 
and  W0  for  each  matrix  stratum.  Another  matrix  must  be  constructed  to  specify  the 
relationship  between  the  values  of  U,  V,  and  W  at  the  upper  and  lower  surfaces  of  a 
matrix  stratum.  The  matrix  shown  below  expresses  U,  V,  and  W  at  the  top  surface  of 
each  matrix  stratum  in  terms  of  U,  V,  and  W  at  the  bottom  of  the  matrix  stratum. 


Expressed  this  way,  the  values  of  U0,  V0,  and  W0  that  are  found  are  equal  to  U,  V,  and 
W  at  the  lower  surface  of  the  matrix  stratum. 
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evaluated  for  matrix  stratum  "1".  The  displacements  for  any  value  of  z  are  found  using 
the  expression 

U(z)  =  e[Q]^  UQ  (51) 

where  is  the  z-coordinate  of  the  lower  surface  of  the  matrix  stratum  in  which  "z"  lies. 

Equation  (49)  is  now  substituted  into  Equation  (46),  so  that  the  external  loads  are 
expressed  in  terms  of  the  unknowns  U0,  V0,  and  W0.  The  coefficients  U0,  V0,  and  W0  are 
found  by  inverting  the  matrix  and  then  multiplying  the  inverse  by  the  force  vector. 
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IV.  A  Simply-Supported  Square  Plate  with  Cross-Plies 


Strata  Theory 

The  strata  theory’s  assumed  displacement  field  does  not,  at  the  outset,  explicitly 
assume  the  x  and  y  dependence  of  u,  v,  and  w.  The  restriction  of  the  problem  to  a 
simply-supported  plate  with  cross-plies  permits  the  use  of  sinusoidal  functions  to  describe 
the  x  and  y  dependencies  in  the  general  PDEs.  The  x  and  y  dependence  can  be  assumed 
to  have  the  form 


w(x,y) 

j4sin(px)sin(£y) 

Vfay) 

Bicos(px)sm(qy) 

C((x,y) 

►  r:  i 

Cisin(px)cos(qy) 

P(x,y) 

9sm(px)sw(qy) 

(52) 


where  p=r/a  and  q=r/b.  The  form  of  the  surface  load  "P"  dictates  this  choice  of  the 
assumed  displacement  field.  These  displacement  field  expressions  and  the  subsequent 
differential  equations  can  be  generalized  for  more  complex  surface  loads  using  Fourier 
analysis. 

Considering  a  composite  with  cross-plies  further  simplifies  the  equations 
previously  developed,  since  it  says  that  S5=S6=S9=0.  When  the  assumed  displacement 
field  expressions  are  substituted  into  Equations  (19),  (20),  and  (21),  the  x  and  y 
dependence  drops  out  entirely.  This  results  in  a  system  of  algebraic  equations  in  the 
unknown  coefficients  A,  B„  and  Q 
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P  =  A  [(jr,0*S,V*<*?*S.°)«2] 

i-0 


(55) 


When  this  is  expressed  in  matrix  form,  this  system  of  equations  can  be  decomposed  into 
a  matrix  of  constants  multiplied  by  the  vector  of  unknown  coefficients  A,  B„  and  Q. 

Using  the  same  displacement  functions  for  u,  v,  and  w  in  Equations  (23)  through 
(28),  the  associated  boundary  conditions  are  now  given  by 
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By  the  definition  of  a  simply  supported  plate,  5u=0  at  y=0  and  y=b,  5v=0  at  x=0  and 
x=a,  and  5w=0  at  all  four  edges.  This  accounts  for  the  boundary  equations  given  in 
Equations  (57),  (59),  (60),  and  (61).  Equations  (56)  and  (58)  are  satisfied  because  <rx  and 
cry  are  defined  to  be  zero  at  the  edges  of  a  simply  supported  plate. 

Specific  Results 

Comparison  calculations  were  made  for  both  isotropic  plates  and  2-ply  composite 
lay-ups.  The  following  values  were  assumed  for  these  test  cases 


Table  1.  Assumed  values  for  comparison  calculations  (square  plate) 


Ef  =  2.0  xlO7 

vf  =  0.35 

volf  =  V2/2 

E*  =  3.0x10s 

=  0.40 

vol»  =  1-V2/2 

p  =  x/a 

a  =  10 

<9  =  1 

q  =  x/b 

b  =  10 

h  =  1 

o 

O 

II 

C 

hn  =  V2/4 

en  =  90° 

hn  =  V2/4 

0.1  =  -(1  +V2)I2 

Pml  =  2W2 

“m2  =  0 

@n2  =  2+V2 

"m3  =  (1  +V'2)/2 

fa  =  2+V2 

Volume  fraction  terms  vol„  and  volf  were  derived  by  assuming  global  volume 
fractions  of  50%  fibers  and  50%  matrix.  Recall  that  vol„  and  volf  refer  to  the  fractions 
of  matrix  and  fiber  in  an  idealized  fiber  ply,  whereas  the  term  "global  volume  fraction" 
refers  to  the  fractions  of  matrix  and  fibers  that  exist  before  the  fiber  plies  are  idealized 
and  dictates  the  thicknesses  of  the  matrix  and  fiber  plies. 
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A  fifth-order  polynomial  in  z  for  the  ST  displacement  field  is  assumed,  i.e.  N=5 
in  Equations  (2)  and  (3),  unless  otherwise  indicated.  Also,  the  applied  load  in  ESP  is 
split  half  and  half  between  the  upper  and  lower  surfaces  of  the  plate.  Loading  the  plate 
this  way  introduces  symmetry  into  the  ESP  solution  that  approximates  the  way  a  load  is 
applied  in  ST. 

The  values  for  h  and  <P  were  chosen  to  normalize  the  governing  equations.  The 
x,  y,  and  z  coordinates,  plate  dimensions  and  displacements  are  normalized  to  the  plate 
thickness  (x/h,  y/h,  z/h,  a/h,  b/h,  u/h,  v/h,  and  w/h).  Also,  the  elastic  moduli  are 
normalized  to  the  applied  pressure  (E^/CP  and  E,/(P). 

The  Isotropic  Square  Plate 

Since  the  plate  consists  of  matrix  only  and  fibers  are  not  present,  Ef  and  vf  are  set 
to  zero  in  the  governing  equations.  When  these  terms  are  set  to  zero  and  N=l,  the  ST 
equations  simplify  to  the  Mindlin  plate  equations  [7:33].  Application  of  the  ESP  approach 
to  an  isotropic  plate  is  the  same  as  for  a  composite,  except  that  no  inter-strata  boundary 
conditions  exists. 

The  complete  set  of  displacements,  stresses  and  strains  are  shown  in  Figures  6. 
through  21.  The  displacements  are  obtained  directly  for  both  ST  and  ESP.  The  strains 
are  derived  from  the  displacement  information,  and  stresses  are  calculated  using  the 
constitutive  stress-strain  relations,  given  in  Equations  (9)  and  (16). 

The  following  plots  for  the  stresses,  strains,  and  displacements  include  only  the 
z-dependent  portion  of  these  parameters.  In  other  words,  the  x-  and  y-dependent  sine  and 
cosine  functions  are  divided  out  for  the  presentation  of  this  data. 
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Figure  9.  ex  for  square  isotropic  plate 
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Figure  10.  ey  for  square  isotropic  plate 
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Figure  11.  cz  for  square  isotropic  plate 
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The  following  tables  list  the  values  of  the  parameters  that  were  calculated  for  the 


preceding  plots: 


Table  2.  Matrix  strata  elastic  coefficients  (square  isotropic  plate) 


Kim 

6.4286E+005 

Kzm 

4.2857E+005 

K3m 

1.0714E+005 

Table  3.  Strata  Theory  coefficients  (square  isotropic  plate) 


A  =  5.3574E-004 

B0  =  0.0 

B,  =  -1.4613E-004 

R2  =  0.0 

B3  =  -2.8843E-005 

B4  =  0.0 

Be  =  -1.7313E-006 


Co  =  0.0 

C,  =  -1.4613E-004 

C2  =  0.0 

C3  =  -2.8843E-005 

C4  =  0.0 

C«  =  -1.7313E-006 


Table  4.  ESP  displacement  function  coefficients  (square  isotropic  plate) 


xu0 

-2.6824E-004 

XV0 

-2.6824E-004 

XW0 

-2.7244E-005 

U0 

-6.6316E-005 

mm 

-6.6316E-005 

W0 

9.0704E-004 

Evaluated  at  z=-h/2 


For  an  isotropic  plate,  there  was  a  reasonable  degree  of  agreement  between  ESP 
and  ST.  There  was  excellent  agreement  for  r„,  ryz,  and  yyt  when  a  cubic  (or  higher 
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order)  polynomial  was  used  for  ST.  The  assumed  order  of  the  polynomial  N  for  ST 
determines  the  variation  of  y «  and  yyz  through  the  thickness.  When  N=l,  ST  predicts 
that  yu,  7yz,  and  ryz  are  constants.  But  as  N  is  increased,  i.e.  N=3  or  larger,  ST  and 
ESP  agree  to  within  1  %  for  these  shear  terms.  The  remaining  stresses  and  strains  as  well 
as  the  displacements  displayed  the  same  general  behavior  for  both  ESP  and  ST,  although 
the  actual  values  calculated  disagreed  by  nearly  a  factor  of  2  in  some  cases. 

The  amount  by  which  ESP  and  ST  disagree  appears  to  be  primarily  dependent 
upon  the  value  of  the  Poisson’s  ratio  (v).  When  the  Poisson’s  ratio  v  is  allowed  to  vary, 
it  is  immediately  apparent  from  Figure  21.  that  ESP  and  ST  agree  exactly  when  v=0, 
and  the  two  solutions  diverge  as  v  increases.  At  i/=0.5,  which  represents  an 
incompressible  material,  ST  predicts  that  w=0,  whereas  ESP  predicts  a  non-zero  z 
displacement.  Intuitively,  the  ESP  results  seem  more  reasonable  when  one  realizes  that 
shear  becomes  the  primary  form  of  deformation  for  an  incompressible  material,  and  some 
displacement  is  still  possible. 

Clearly,  the  effect  of  shear  is  diminished  in  ST.  ST  assumes  that  ez=0,  while  ESP 
expects  a  nonzero  ez.  Since  the  normal  stresses  are  coupled  through  v,  elimination  of  one 
of  the  normal  stresses,  in  this  case  ez,  introduces  an  error  that  is  not  present  in  ESP. 

The  2-Ply  Square  Composite  Plate 

The  complete  set  of  displacements,  stresses  and  strains  are  shown  in  Figures  22. 
through  36.  The  displacements  are  obtained  directly  from  both  ST  and  ESP.  The  strains 
are  derived  from  the  displacement  information,  and  stresses  are  calculated  using  the 
constitutive  stress-strain  relations,  given  in  Equations  (9)  and  (16). 
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The  following  plots  for  the  stresses,  strains,  and  displacements  include  only  the 
z-dependent  portion  of  these  parameters.  In  other  words,  the  x-  and  y-dependent  sine  and 
cosine  functions  are  divided  out  for  the  presentation  of  this  data. 


-e-  U  -  E.S.P.  -**-  U  -  Strata  Theory 

Figure  22.  In-plane  displacement  u  for  square  2-ply  composite  plate 
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-e-  V  -  E.S.P.  -**—  V  -  Strata  Theory 

Figure  23.  In-plane  displacement  v  for  square  2-ply  composite  plate 


s-  E.S.P.  x  Strata  Theory 

Figure  24.  Out-of-plane  displacement  w  for  square  2-ply  composite  plate 
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Normalized  Z  Location 


Strain  -  Gamma  XY 
(Times  10E-8) 

-e-  E.S.P.  Strata  Theory 

Figure  28.  y v  for  square  2-ply  composite  plate 
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Normalized  Z  Location 


0.5< 


(Times  10E-6) 


-e-  E.S.P.  —**—  Strata  Theory 

Figure  29.  ya  for  square  2-ply  composite  plate 


-e-  ES.P.  — Strata  Theory 

Figure  30.  yn  for  square  2-ply  composite  plate 
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E.S.P.  -**-  Strata  Theory 

Figure  36.  rn  for  square  2-ply  composite  plate 


52 


The  following  tables  list  the  values  of  the  parameters  that  were  calculated  for  the 


preceding  plots: 


Table  5.  Matrix  strata  elastic  coefficients  (2-ply  square  plate) 


Klm 

6.4286E+005 

K2m 

4.2857E+005 

1.0714E+005 

Table  6.  Fiber  strata  elastic  coefficients  (2-ply  square  plate) 


Fiber  Stratum  1  Fiber  Stratum  1 


Si, 

2.2886E+007 

2.1949E+006 

S2f 

2.1949E+006 

2.2886E+007 

| 

1.2347E+007 

1.2347E+007 

$4f 

3.6581E+005 

3.6581E+005 

Sjf 

O.OOOOE+OOO 

O.OOOOE+OOO 

mm 

O.OOOOE+OOO 

O.OOOOE+OOO 

S7f 

5.2692E+006 

3.6581E+005  1 

S.f 

3.6581E+005 

5.2692E+006 

mm 

O.OOOOE+OOO 

O.OOOOE+OOO 

Siof 

1.2347E+007 

4.2857E+005 

Suf 

4.2857E+005 

1.2347E+007 

Si2f 

2.2886E+007 

2.2886E+007 

I  Sur 

O.OOOOE+OOO 

O.OOOOE+OOO 

When  fibers  are  present,  the  agreement  between  ST  and  ESP  for  the 
displacements  u,  v,  and  w  are  very  good.  Also,  there  is  very  good  agreement  between 


Table  7.  Strata  Theory  coefficients  (2-ply  square  plate) 


A  =  2.5179E-006 
B0  =  -3.1637E-007 
B,  =  1.0810E-007 
B2  =  -3.2442E-010 
Bj  =  -1.3121E-006 
B4  =  2. 1974E-009 
Bs  =  2.1453E-006 


C0  =  3.1637E-007 
C,  =  1.0810E-007 
C2  =  3.2442E-010 
C3  =  -1.3121E-006 
C4  =  -2. 1974E-009 
C,  =  2.1453E-006 


Table  8.  ESP  displacement  function  coefficients  (2-ply  square  plate) 


Matrix  Stratum  1  Matrix  Stratum  2  Matrix  Stratum  3 


XUo 

-7. 1529E-007 

1.0012E-007 

-8.0927E-007 

AV0 

-7.2279E-007 

9.2623E-008 

-8.0927E-007 

AW0 

7.8046E-007 

-9.7706E-009 

-7.573  IE-007 

U0 

-2.5245E-007 

-2.6659E-007 

-2.1051E-007 

n 

2.6659E-007 

2.5245E-007 

3.0825E-007 

Wo 

2.5197E-006 

2.5197E-006 

2.5760E-006 

Evaluated  at  the  z  coordinate  at  the  bottom  of  each  matrix  stratum 


ESP  and  ST  in  the  general  behavior  for  all  the  stresses  and  strains,  with  the  obvious 
exception  of  ez.  The  polynomial  form  of  the  z  dependence  of  the  assumed  displacement 
field  is  clearly  evident  in  the  plots  of  the  shear  stresses  and  shear  strains. 

In  comparing  the  results  of  the  2-ply  composite  with  the  isotropic  plate,  it  is 
obvious  that  the  fibers  dominate  the  behavior  of  the  composite.  While  the  general 
behavior  of  the  ST  results  echo  those  of  ESP,  there  is  still  some  disagreement  in  the 
actual  values  calculated  for  displacements,  stresses,  and  strains.  Other  test  cases  for  a 
square  plate  were  examined,  and  they  revealed  that  the  level  of  agreement  between  ESP 
and  ST  seems  to  be  most  highly  dependent  upon  the  Poisson’s  ratio  for  the  fibers  (i>f), 
and  the  order  of  the  polynomial  (N)  assumed  for  ST. 
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Just  as  with  the  isotropic  plate,  there  is  considerable  disagreement  in  ez  between 
ST  and  ESP  for  a  composite  panel,  but  this  should  be  expected.  In  the  matrix  strata  that 
comprise  the  top  and  bottom  surfaces  of  the  composite,  note  that  ez  increases  slightly  as 
z  gets  farther  away  from  the  free  surfaces  at  z=±h/2.  This  demonstrates  a  coupling  of 
the  normal  stresses.  At  the  free  surfaces  of  the  composite,  the  matrix  is  allowed  to 
deform  freely  through  Poisson  effects.  However,  at  the  matrix-fiber  strata  interface,  the 
matrix  is  forced  to  match  the  displacements  of  the  fibers,  and  so  it  cannot  deform  freely. 

Because  ez  is  allowed  to  be  nonzero  for  ESP,  ESP  predicts  that  the  presence  of 
fibers  results  in  significant  jumps  in  ez  through  the  thickness.  Through-the-thickness  shear 
through  the  fiber  strata  affects  how  ez  varies  for  the  ESP  solution.  Since  a  square  plate 
is  being  considered  here,  and  since  the  plies  are  rotated  90°  with  respect  to  each  other, 
dimensional  symmetry  causes  the  fiber  strata  to  absorb  equal  amounts  of  strain  energy. 
If  the  aspect  ratio  of  the  plate  (a/b)  is  allowed  to  increase,  more  strain  energy  is  absorbed 
by  the  fiber  stratum  whose  fibers  are  aligned  with  the  shorter  dimension  of  the  plate  than 
by  a  stratum  whose  fibers  are  parallel  with  the  long  dimension  of  the  plate.  This  will  be 
more  clearly  evident  in  the  cylindrical  bending  case  that  follows. 

Assuming  that  the  aspect  ratio  b/a  is  constant,  both  ST  and  ESP  predict  that  the 
z  displacement  (w)  of  a  composite  varies  as  a  function  of  the  plate  width  squared  (a2). 
In  other  words,  the  limit  of  w/a2  as  "a"  goes  to  infinity  is  a  constant.  When  the  given 
load  distribution  is  used,  and  a  linear  polynomial  is  assumed  for  ST  (i.e.  N=l),  the 
solution  for  the  coefficient  "A"  of  the  ST  assumed  displacement  field  predicts  this  limit 
to  be 
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lim  —  = 

«-•  a 2 


p 


*2[s°f +  sl\ 


(62) 


This  is  a  fundamentally  different  behavior  than  that  of  an  isotropic  plate  or  for  a  set  of 
stacked  orthotropic  plates.  CPT,  the  exact  orthotropic  plate  solutions  developed  by 
Pagano  [8:398-411;  9:20-34],  ESP,  or  ST  (when  applied  to  an  isotropic  plate)  predicts 
that  the  limit  of  w/a4  as  "a"  goes  to  infinity  is  a  constant.  Note  the  dependence  on  the 
characteristic  plate  width  "a".  Using  the  Kirchhoff  plate  assumptions  for  a  thin  isotropic 
plate  (Navier’s  Solution),  this  limit  is  given  as  [11:548-549] 


Dm  *  = 

«—  a 4  /j3tc4£ 


(63) 


When  "a"  is  allowed  to  get  large,  the  ESP  results  corresponds  to  this 
the  material  property  terms  correspond  to  what  would  be  a  K,  term 
elasticity  problem. 

The  expression  for  an  isotropic  plate  using  ST  becomes 

jjjh  A_  __  3P(l-2v„Xl+v.) 

a 4  h*n*Em(l-vm) 


value.  Notice  that 
for  a  plane  stress 


(64) 


Similar  expressions  arise  for  CPT  or  the  exact  orthotropic  plate  solutions  as  well.  Notice 
that  the  material  property  terms  corresponds  to  an  expression  for  a  Kt  term  that  would 
be  used  for  a  plane  strain  elasticity  problem. 
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V.  The  Cylindrical  Bending  Case 


Cylindrical  bending  is  a  special  case  of  the  rectangular  plate  formulation 
developed  above.  The  dimension  in  the  y  direction  (given  as  "b")  is  assumed  to  be 
infinite,  therefore  v  and  derivatives  with  respect  to  y  are  zero.  As  a  result,  some  strain 
terms  are  defined  as  being  zero.  Because  "b"  is  infinite,  ey,  r^,  and  ryz  are  zero 
throughout  the  composite. 

Strata  Theory 

Using  the  fact  that  b  =  oo  t  and  that  v  and  derivatives  with  respect  to  y  equal  0,  the 
PDEs  in  Equations  (19)  and  (21)  reduce  to 

ft.**'*,  -  E  { }  j- wu(  (6S) 

i*0 

and 

(£° +$7)w;a+5^  =0  (66) 

i-0 


where 


m  Zml 

*,'=E  /(«.,* 


/-I 


*»/ 


*-l 


(67) 


.7 


These  are  the  same  definitions  as  used  for  the  rectangular  plate.  Note  that  Equation  (20) 
vanishes  for  the  cylindrical  bending  case. 

The  boundary  conditions  that  remain  are 

N 

E 

t-0 


M  *■<  JV 


Zml 


5w 


=  0  <«» 


Simply-Supported  Cylindrical  Bending  with  Cross-Plies 

For  simply-supported  cylindrical  bending,  the  x  dependence  of  w,  %,  and  P  can 
be  assumed  to  have  the  following  form 


w(x) 

A  sin (px) 

%(x) 

►  =  - 

B{cos(px) 

Pix ) 

Psin  (px) 

(70) 


where  p=x/a. 

These  expressions  are  now  used  in  the  above  PDEs,  which  results  in  a  system  of 
equations  in  the  unknown  coefficients  A  and  B; 

o  =  AutmKi  'P]  ♦£ «,  [  wIk‘;‘-k  pY;‘*k['’)  ]  m> 

r  0,1,2,  A 
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and 


P  =  A(K2+S°)p2  +  J2  BlWmK'i'p]  (72) 


When  this  is  expressed  in  matrix  form,  this  system  of  equations  can  be  decomposed  into 
a  matrix  of  constants  multiplied  by  the  vector  of  unknown  coefficients  A  and  B;. 


Exact  Stratified  Plate 

As  with  the  ST  approach,  cylindrical  bending  is  a  special  case  of  the  rectangular 
plate  formulation  developed  above.  Again,  the  dimension  in  the  y  direction  (given  as  "b") 
is  assumed  to  be  infinite,  therefore  v  and  derivatives  with  respect  to  y  are  zero.  The 
assumed  displacement  field  expressions  are 


uM  | 

j  U(z)cos(px) 

wfoz)  j 

\  Rr(z)sin(px) 

(73) 


where  p=r/a. 

When  v  and  the  y  derivatives  are  set  to  zero  in  the  rectangular  plate  formulation, 
and  b  is  defined  to  be  infinite,  the  ODEs  that  result  are 


*3.  0 

\u"\  + 

o  pKAm 

M  + 

o 

■ 

If 

fl/1 

fol 

0  K . 

W 

-pK4m  0 

w 

0  -p% 

*  J 

Vj 

10 1 

where  KH=K2-t-K3.  When  the  2  second-order  ODEs  are  recast  as  a  set  of  four  first-order 
coupled  ODEs,  U  and  W  are  assumed  to  have  the  form  of  Uge^,  and  the  governing 
equations  are  expressed  as  a  system  of  4  coupled  algebraic  equations.  Rewritten  in  state- 
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space  form,  they  become 


where 


[K?]*M*]]{ffo}-<o> 


{*.}- 


p% 


[<?]  -  -- 


p%. 


The  eigenvalues  of  this  system  are  2  sets  of  repeated  roots 


k  =  ±  p 


The  same  method  for  generating  an  eigenvector  matrix  as  described  above  for  the 
rectangular  plate  is  used  [1:53-57],  but  it  is  tailored  for  a  system  with  four  eigenvalues. 
First,  the  functions  f(z,X)  and  g(z,X)  are  defined  as 

g(z,X)  =  o0(z)  +  a  j(z)  k  +  a2(z)  k2  +  a3(z)  X3 
f(z,k)  =  exz 


and  these  expressions  are  manipulated  to  solve  for  the  a;’s  as  before. 


1  X  X2 

X3 

f(z,X) 

e * 

3X2 

ao 

Z'feA) 

►  —  < 

ze ** 

►  =: 

0  1  2X 

i 

/fc-« 

e'Xl 

1  -X  X2 

-X3 

«2 

ze  ~kz 

3X2 

o 

(80) 


Finally,  exponentiation  of  the  matrix  [Q]  can  be  obtained  by  the  relation 

*1<?k  -  «0[/l  +  *JLQ] +  «2[<?]2  +  «3[<?]3  <81> 

The  boundary  conditions  are  constructed  in  the  same  fashion  as  for  the  rectangular 
plate,  using  a  series  of  force  balance  equations  for  the  x-z  shear  and  normal  stress. 
Equilibrium  is  established  at  die  upper  surface,  at  each  fiber  stratum,  and  finally  at  the 
lower  surface. 


The  boundary  conditions  are  setup  the 

same  as 

for  the  rectangular  plate,  given  by 

Mlx 

M2, 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-Mi, 

-M2, 

Mi2 

M2n 

0 

0 

0 

0 

0 

0 

0 

0 

0 

/ 

0 

-/ 

0 

0 

0 

0 

0 

0 

[BC]  = 

0 

0 

0 

0 

0 

0 

-Mi2 

-M22 

Mi3 

M2n 

0 

0 

0 

0 

0 

0 

0 

0 

0 

I 

0 

-I 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-Mi3 

-M23 

(82) 

but  the  submatrices 

are  defined  as 
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Ml,= 


K .  0 

imi 

0  K , 


M2,= 


M2/k  = 


V  A 

'.1  \ 

0 

pK,m 

0 

p\ 

0 

0  P2S7<1 


V  +  ^k+i 


instead.  These  represent  stresses  r„  and  <r2  at  the  top  and  bottom  surfaces  of  the 
composite  stack,  and  at  the  fiber  strata. 


Specific  Results 

Comparison  calculations  were  made  for  both  isotropic  plates  and  2-ply  composite 
lay-ups.  The  following  values  were  assumed  for  these  test  cases 
Table  9.  Assumed  values  for  comparison  calculations  (cylindrical  bending) 


Ef  =  10* 


Em  =  107 


p  =  r/a 


=  0 


9n  =  0° 


ef2  =  90° 


=  <l+V2)/2 


(*02  —  0 


a„3  =  (lW2)/2 


vf  =  0.25 


a  =  10 


00 


volf  =  V'2/2 


vol^  =  1V2/2 


<P  =  1 


=  1 


0mi  =  2W2 


0*2  =  2W2 


0*3  =  2W2 
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Volume  fraction  terms  voln,  and  volf  were  derived  by  assuming  global  volume 
fractions  of  50%  fibers  and  50%  matrix.  Recall  that  vol,,,  and  volf  refer  to  the  fractions 
of  matrix  and  fiber  in  an  idealized  fiber  ply,  whereas  the  term  "global  volume  fraction" 
refers  to  the  fractions  of  matrix  and  fibers  that  exist  before  the  fiber  plies  are  idealized 
and  dictates  the  thicknesses  of  the  matrix  and  fiber  plies. 

A  fifth-order  polynomial  in  z  for  the  ST  displacement  field  is  assumed,  i.e.  N=5 
in  Equations  (2),  unless  otherwise  indicated.  Also,  the  applied  load  is  split  half  and  half 
between  the  upper  and  lower  surfaces  of  the  plate. 

The  values  for  h  and  <P  were  chosen  to  normalize  the  governing  equations.  The 
x  and  z  coordinates  and  the  plate  width  are  normalized  to  the  plate  thickness  (x/h,  z/h, 
and  a/h),  and  the  elastic  moduli  are  normalized  to  the  applied  pressure  (E^/tP  and  E,/(P). 
Displacements  are  also  normalized  to  the  plate  thickness  (u/h  and  w/h). 

The  Isotropic  Cylindrical  Bending  Case 

Since  the  plate  consists  of  matrix  only  and  fibers  are  not  present,  Ef  and  vf  are  set 
to  zero  in  the  governing  equations.  The  complete  set  of  displacements,  stresses  and 
strains  are  shown  in  Figures  37.  through  45.  The  displacements  are  obtained  directly 
from  both  ST  and  ESP.  The  strains  are  derived  from  the  displacement  information,  and 
stresses  are  calculated  using  the  constitutive  stress-strain  relations,  given  in  Equations  (9) 
and  (16). 

The  following  plots  for  the  stresses,  strains,  and  displacements  include  only  the 
z-dependent  portion  of  these  parameters.  In  other  words,  the  x-  and  y-dependent  sine  and 
cosine  functions  are  divided  out  for  the  presentation  of  this  data. 


63 


-0.00020  -0.00015  -0.00010  -0.00005  -0.00000  0.00005  0.00010  0.00015  0.00020 
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-B-  E.S.P.  -**-  Strata  Theory 

Figure  41.  ya  for  cylindrical  oending  isotropic  plate 


o  E.S.P.  "X  Strata  Theory 

Figure  42.  <rz  for  cylindrical  bending  isotropic  plate 
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-®~  E.S.P.  -**-  Strata  Theory 

Figure  45.  Displacement  w  for  cylindrical  bending  isotropic  plate  vs.  v 


The  following  tables  list  the  values  of  the  parameters  that  were  calculated  for  the 
preceding  plots: 


Table  10.  Matrix  strata  elastic  coefficients  (cylindrical  bending  isotropic  plate) 


K,. 

6.4286E+005 

K* 

4.2857E+005 

K- 

1.0714E+005 

Just  as  with  the  rectangular  plate,  there  was  a  reasonable  degree  of  agreement 
between  ESP  and  ST.  There  was  excellent  agreement  for  r„,  ax,  and  yu  when  a  cubic 
(or  higher  order)  polynomial  was  used  for  ST.  The  remaining  stresses  and  strains  as  well 
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Table  11.  Strata  Theory  coefficients  (cylindrical  bending  isotropic  plate) 

A  =  2.0297E-003 
B0  =  0.0 

=  -5.9320E-004 

Bj  =  0.0 

Bj  =  -5.8545E-005 

B4  =  0.0 

Bs  =  -1.7452E-006 

Table  12.  ESP  displacement  function  coefficients  (cylindrical  bending  isotropic  plate) 


XU0 

-1.1051E-003 

xw0 

1.1245E-004 

U0 

5.4064E-004 

W0 

3.5178E-003 

Evaluated  at  z=-h/2 


as  the  displacements  displayed  the  same  general  behavior  for  both  ESP  and  ST,  although 
the  actual  values  calculated  disagreed  by  nearly  a  factor  of  2  in  some  cases. 

The  amount  by  which  ESP  and  ST  disagree  appears  to  be  primarily  dependent 
upon  the  value  of  the  Poisson’s  ratio  (r).  When  the  Poisson’s  ratio  v  is  allowed  to  vary, 
it  is  immediately  apparent  from  Figure  45.  that  ESP  and  ST  agree  exactly  when  v=0, 
and  the  two  solutions  diverge  as  v  increases. 

Just  as  with  the  rectangular  plate,  ST  assumes  that  ez=0,  while  ESP  expects  a 
nonzero  ez.  Since  the  normal  stresses  are  coupled  through  v,  elimination  of  ez  introduces 
an  error  that  is  not  present  in  ESP.  Again,  and  there  is  poor  agreement  for  fz  and  <rz. 
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The  2-Ply  Cylindrical  Bending  Case 

The  complete  set  of  displacements,  stresses  and  strains  are  shown  in  Figures  46. 
through  53.  The  displacements  are  obtained  directly  from  both  ST  and  ESP.  The  strains 
are  derived  from  the  displacement  information,  and  stresses  are  calculated  using  the 
constitutive  stress-strain  relations,  given  in  Equations  (9)  and  (16). 

The  following  plots  for  the  stresses,  strains,  and  displacements  include  only  the 
z-dependent  portion  of  these  parameters.  In  other  words,  the  x-  and  y-dependent  sine  and 
cosine  functions  are  divided  out  for  the  presentation  of  this  data. 


s-  E.S.P.  -**-  Strata  Theory 

Figure  46.  In-plane  displacement  u  for  cylindrical  bending  2-ply  composite  plate 
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Displacement  Function  W(z) 
(Times  10E-6) 


-e-  E.S.P.  -**-  Strata  Theory 

Figure  47.  Out-of-plane  displacement  w  for  cylindrical  bending  2-ply  composite  plate 


-e-  E.S.P.  — Strata  Theory 

Figure  48.  et  for  cylindrical  bending  2-Ply  composite  plate 


Normalized  Z  Location 


Stress  -  Sigma  X 


-e~  E.S.P.  Strata  Theory 

Figure  51.  ax  for  cylindrical  bending  2-Ply  composite  plate 


-e-  E.S.P.  -**-  Strata  Theory 

Figure  52.  at  for  cylindrical  bending  2-Ply  composite  plate 


-e-  E.S.P.  -**-  Strata  Theory 

Figure  53.  ru  for  cylindrical  bending  2-Ply  composite  plate 


The  following  tables  list  the  values  of  the  parameters  that  were  calculated  for  the 
preceding  plots: 


Table  13.  Matrix  strata  elastic  coefficients  (2-ply  cylindrical  bending  plate) 


6.4286E+005  I 

K2m 

4.2857E+005 

8 

1.0714E+005 

When  fibers  are  present,  the  agreement  between  ST  and  ESP  for  the 
displacements  u,  v,  and  w  are  very  good.  Also,  there  is  very  good  agreement  between 
ESP  and  ST  in  the  general  behavior  for  all  the  stresses  and  strains,  with  the  obvious 
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Table  14.  Fiber  strata  elastic  coefficients  (2-ply  cylindrical  bending  plate) 

Fiber  Stratum  1  Fiber  Stratum  1 


S,f 

2.2886E+007 

2.1949E+006 

mm 

2.1949E+006 

2.2886E+007 

I 

1.2347E+007 

1.2347E+007 

□z 

3.6581E+005 

3.6581E+005 

mm 

O.OOOOE+OOO 

O.OOOOE+OOO 

O.OOOOE+OOO 

O.OOOOE+OOO 

1 _ s?f 

5.2692E+006 

3.6581E+005 

WM 

3.6581E+005 

5.2692E+006 

O.OOOOE+OOO 

O.OOOOE+OOO 

Si« 

1.2347E+007 

4.2857E+005 

S„f 

4.2857E+005 

1.2347E+007 

S|2f 

2.2886E+007 

2.2886E+007 

Sl3f 

O.OOOOE-fOOO 

O.OOOOE+OOO 

Table  15.  Strata  Theory  coefficients  (2-ply  cylindrical  bending  plate) 

A  =  5.0817E-006 
B0  =  8.6071E-008 
B,  =  -4. 1525E-007 
Bj  =  2.6488E-010 
B3  =  -1.9537E-007 
B4  =  -1.7935E-009 
Bs  =  3.1359E-007 

exception  of  ez.  The  polynomial  form  of  the  z  dependence  of  the  assumed  displacement 
field  is  clearly  evident  in  the  plots  of  the  shear  stresses  and  shear  strains. 

In  comparing  the  results  of  the  2-ply  composite  with  the  isotropic  plate,  it  is 
obvious  that  the  fibers  dominate  the  behavior  of  the  composite.  While  the  general 
behavior  of  the  ST  results  echo  those  of  ESP,  there  is  still  some  disagreement  in  the 
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Table  16.  ESP  displacement  function  coefficients  (2-ply  cylindrical  bending  plate) 


Matrix  Stratum  1  Matrix  Stratum  2  Matrix  Stratum  3 


XU0 

-1.5197E-006 

-1.4293E-006 

-1.6406E-006 

xw0 

7.7067E-007 

-6.4056E-007 

-7.1890E-007 

U0 

-3.3253E-008 

1.6444E-007 

2.8112E-007 

W0 

5.0720E-006 

5.1687E-006 

5.2223E-006 

Evaluated  at  the  z  coordinate  at  the  bottom  of  each  matrix  stratum 


actual  values  calculated  for  displacements,  stresses,  and  strains.  Other  test  cases  for  a 
square  plate  were  examined,  and  they  revealed  that  the  level  of  agreement  between  ESP 
and  ST  seems  to  be  most  highly  dependent  upon  the  Poisson’s  ratio  for  the  fibers  (i>f), 
and  the  order  of  the  polynomial  (N)  assumed  for  ST. 

Just  as  with  the  isotropic  plate,  there  is  considerable  disagreement  in  ez  between 
ST  and  ESP  for  a  composite  panel,  but  this  should  be  expected.  In  the  matrix  strata  that 
comprise  the  top  and  bottom  surfaces  of  the  composite,  note  that  ez  increases  slightly  as 
z  gets  farther  away  from  the  free  surfaces  at  z=±h/2.  This  demonstrates  a  coupling  of 
normal  stresses.  At  the  free  surfaces  of  the  composite,  the  matrix  is  allowed  to  deform 
freely  through  Poisson  effects.  However,  at  the  matrix-fiber  strata  interface,  the  matrix 
is  forced  to  match  the  displacements  of  the  fibers,  and  so  it  cannot  deform  freely. 

Because  ez  is  allowed  to  be  nonzero  for  ESP,  ESP  predicts  that  the  presence  of 
fibers  results  in  significant  jumps  in  ez  through  the  thickness.  As  indicated  for  the 
rectangular  plate  case,  the  through-the-thickness  shear  through  the  fiber  strata  affects  how 
ez  varies  for  the  ESP  solution.  A  fiber  ply  where  the  fibers  span  the  width  of  the  plate 
will  absorb  much  more  strain  energy  than  a  ply  with  fibers  that  run  parallel  to  the 
simply-supported  edges.  This  induces  a  larger  jump  in  €z  between  the  top  and  middle 
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matrix  strata  than  between  the  middle  and  bottom  matrix  strata  [Figure  49.].  If  the  fiber 
plies  were  reversed  (90  &  0  instead  of  0  &  90),  there  would  be  a  greater  difference  in 
ez  across  the  second  fiber  ply  than  the  first,  and  ez  through  the  middle  matrix  stratum 
would  be  positive  by  the  same  average  magnitude  (i.e.  it  would  be  "flipped"  to  the 
positive  side  of  the  graph).  The  cylindrical  bending  case  is  the  limiting  case  of  a 
rectangular  plate  with  a  high  aspect  ratio. 

As  was  indicated  for  a  rectangular  plate,  both  ST  and  ESP  predict  that  the  z 
displacement  (w)  of  a  composite  varies  as  a  function  of  the  plate  width  squared  (a2).  In 
other  words,  the  limit  of  w/a2  as  "a"  goes  to  infinity  is  a  constant.  When  the  given  load 
distribution  is  used,  and  a  linear  polynomial  is  assumed  for  ST  (i.e.  N=l),  the  solution 
for  the  coefficient  "A"  of  the  ST  assumed  displacement  field  predicts  this  limit  to  be 


lim 

a-» 


p 


(84) 


This  is  a  fundamentally  different  behavior  than  that  of  an  isotropic  plate  or  for  a  set  of 
stacked  orthotropic  plates.  CPT,  the  exact  orthotropic  plate  solutions  developed  by 
Pagano  [8:398-411;  9:20-34],  ESP,  or  ST  (when  applied  to  an  isotropic  plate)  predicts 
that  the  limit  of  w/a4  as  "a"  goes  to  infinity  is  a  constant.  Note  the  dependence  on  the 
characteristic  plate  width  "a".  Applying  Navier’s  solution  to  the  cylindrical  bending  case 
for  an  isotropic  material,  the  normalized  displacement  is 


12P(l-v2) 


w 


(85) 
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The  expression  for  an  isotropic  plate  using  ST  becomes 


a-m 


a 


4 


12P(l-2vw)(l^vw) 

A*«4£.(1-v.) 


(86) 


Similar  expressions  arise  for  CPT  or  the  exact  orthotropic  plate  solutions  as  well. 
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V7.  Conclusions  and  Recommendations 


Interpretation  of  Results 

For  an  isotropic  plate,  there  was  a  reasonable  degree  of  agreement  between  ESP 
and  ST.  There  was  excellent  agreement  for  ru,  ryz,  ya,  and  yyl  when  a  cubic  (or  higher 
order)  polynomial  was  used  for  ST.  The  remaining  stresses  and  strains  as  well  as  the 
displacements  displayed  the  same  general  behavior  for  both  ESP  and  ST,  although  the 
actual  values  calculated  disagreed  by  nearly  a  factor  of  2  in  some  cases.  Since  the  normal 
stresses  are  coupled  through  p,  and  since  ez  is  set  to  zero  in  ST,  this  introduces  an  error 
that  is  not  present  in  ESP.  As  a  result,  the  amount  by  which  ESP  and  ST  disagree  is 
heavily  dependent  on  the  value  of  the  Poisson’s  ratio  (r). 

When  fibers  are  present,  however,  the  agreement  between  ST  and  ESP  for  the 
displacements  u,  v,  and  w  are  very  good.  Also,  there  is  very  good  agreement  between 
ESP  and  ST  in  the  general  behavior  for  all  the  stresses  and  strains,  with  the  obvious 
exception  of  et.  While  it  is  clear  that  the  fibers  dominate  the  behavior  of  the  composite, 
the  level  of  agreement  in  the  actual  values  calculated  seems  to  be  most  highly  dependent 
upon  the  Poisson’s  ratio  for  the  fibers  (r{). 

Since  the  difference  in  the  displacements  between  the  ESP  and  ST  solutions  for 
a  composite  are  very  small,  the  approximation  that  epsilon  z  is  zero  appears  to  be  a 
reasonably  valid  assumption  for  ST.  Although  ST  was  not  really  intended  to  be  used  for 
detailed  stress/strain  prediction,  the  fact  that  ST  and  ESP  agree  so  well  is  another 
confirmation  of  the  general  validity  of  the  ST  approach.  It  also  reinforces  the  potential 
of  ST  as  a  dynamic  analysis  tool. 
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The  applied  load  used  for  the  example  calculations  had  a  simple  sinusoidal 
distribution,  but  Fourier  analysis  can  be  used  to  decompose  any  arbitrary  load  into  an 
infinite  sum  of  sines  and  cosines.  The  equations  used  for  the  example  problems  can 
readily  be  generalized  for  more  complex  loads  by  modifying  the  values  of  p  and  q. 
However,  if  many  terms  of  the  Fourier  series  are  needed  to  represent  the  applied  load, 
considerable  computation  time  may  be  required. 

Areas  for  Further  Research 

It  is  not  clear  whether  the  behavior  of  the  normalized  displacement  w  as  a 
function  of  plate  width  represents  the  true  behavior  of  a  composite  panel.  This  would 
seem  to  be  an  important  characteristic  of  the  Stratified  Plate  Model  that  needs  to  be 
verified,  or  refuted,  experimentally. 

The  ST  equation  development  readily  lends  itself  to  a  formulation  that  includes 
a  dynamic  term.  If  the  ST  method  is  to  be  used  for  mode  shape  prediction,  verification 
of  the  behavior  of  ST  with  respect  to  the  plate  dimensions,  as  suggested  above,  is  needed 
as  a  preliminary  and  preparatory  step  to  determine  the  suitability  of  ST  as  a  vibration 
analysis  tool. 

There  are  indications  that  the  ST  formulation  would  lend  itself  well  to  the 
development  of  a  new  type  of  finite  element.  Indeed,  some  preliminary  work  in  this  area 
has  already  been  done.  Once  the  formulation  of  such  an  element  is  developed,  validation 
and  verification  is  needed  to  determine  such  aspects  as  convergence  and  compatibility 
within  and  between  elements. 

One  of  the  strengths  of  the  ST  method  is  that  the  order  of  the  polynomial  (N)  for 
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the  assumed  z  dependence  in  ST  can  be  selected  as  necessary  to  obtain  the  desired  level 
of  accuracy.  This  study  did  not  explicitly  examine  the  effect  of  N  on  the  accuracy  of  the 
results,  particularly  for  multi-ply  composites,  but  it  is  likely  that  the  value  of  N  that 
yields  the  best  trade-off  between  computation  time  and  accuracy  will  be  a  matter  of 
experience  and  the  specific  details  of  the  problem  being  studied.  Nevertheless,  ST  should 
be  applied  to  multi-ply  composites  to  ensure  convergence  and  determine  the  order  of  the 
polynomial  (N)  required  for  adequate  accuracy. 

Also,  ST  should  be  applied  to  lay-ups  with  angle  plies  to  examine  the  coupling 
effects,  and  to  determine  possible  forms  of  solutions  for  such  cases.  The  x-  and  y- 
dependent  portions  of  the  assumed  displacement  field  for  a  simply-supported  rectangular 
plate  with  cross-plies  were  fairly  easy  to  obtain,  but  the  form  of  the  displacement  field 
expressions  when  angle-plies  are  present  are  not  so  obvious  and  should  be  explored 


further. 


Appendix  A:  MATLAB  Routines  for  Finite  Rectangular  Composite  Panel 


SOLVE 

A=10 

B*10 

setparam 

calparam 

exactsol 

■txplies 

quit 


SETPARAM 

%  SETPARAM. M  routine  for  PC-MATLAB 

% 

format  short  e 

global  n  m  A  B  p  q  matprops  fibprops  zdist  N  alpha  beta  tf  tm  FM  R  RMI 
diary  setparam.log 
% 

iso*0 

n»l 

m»l 

« 

Em*0 . 3e6 
NUm* . 40 

Gm-Em/(l+NUm)/2 

% 

kl-Em* ( 1-NUm) / ( l-2*NUm) / ( l+WUm) 
k2-Em*NUm/ (l-2*NUm) / ( 1+NUm) 
k3*Gm 

% 

KM* I kl  k2  k2  0  0  0; 

k2  kl  k2  0  0  0; 

k2  k2  kl  0  0  0; 

0  0  0  k3  0  0; 

0  0  0  0  k3  0; 

0  0  0  0  0  k3  ] ; 

% 

matprops* [kl  k2  k3] 

% 

if  iso«*l 

tf * [ 0  0] 
tm*(.25  .5  .25] 
alpha*0 
beta*l 
« 

fibprops* [0  00000000; 

00C000000] 

else 

sqrt2=sqrt(2) 

volf*sqrt2/2 

volm*l-volf 

% 

tf*[sqrt2/4  sqrt2/4] 

tm*[ .25-sqrt2/8  ,5-sqrt2/4  .25-sqrt2/8) 

alpha*- ( l+sqrt2 ) /2 

beta*2+sqrt2 

« 
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Ef=20e6 


NUf = .  35 

Gf -Ef / ( 1+NUf ) / 2 

% 

al-volf *Ef * ( 1-NUf ) / ( l-2*NUf ) / ( 1+NUf ) +volm*kl 
B2=kl/volm 

s3-volf *Ef *NUf / ( l-2*NUf ) / ( 1+NUf ) +volm*k2 

84- Gm/volm 

85- 0 

86- 0 

b7=vo 1 f  *Gf +vo lm*Gm 

s8-Gm/volm 

a9-0 

sl0-k2 


% 

fibprop8=(al  s2  s3  b4  s5  b6 
s2  si  s3  s4  s5  s6 


« 

SMI- [si  83 
s3  b2 
b3  810 
0  0 
0  0 
0  0 

« 

SM2-[s2  b3 
b3  b1 
BlO  83 
0  0 
0  0 
0  0 

end 


0  0  0; 

0  0  0; 

0  0  0; 

a4/2  0  0; 

0  87/2  0; 

0  0  s8/2  ] 

0  0  0; 

0  0  0; 

0  0  0; 

84/2  0  0; 

0  s8/2  0; 

0  G  87/2] 


« 

zdit  t- [ . 25 ,  -.25] 
N-5 


b7  b8  b9 ; 
b8  87  s9 ] 


CALPARAM 

%  CALPARAM. M  routine  for  PC-MATLAB 

% 

p-n*pi/A 

q*m*pi/B 

% 

R2-p*p+q*q 
R-sqrt (R2 ) 

% 


RM*  [  1 

R 

R2 

R*R2 

R2+R2 

R*R2*R2; 

0 

1 

2*R 

3*R2 

4*R*R2 

5*R2*R2 ; 

0 

0 

2 

6*R 

12+R2 

20*R*R2; 

1 

-R 

R2 

-R*R2 

R2*R2 

-R*R2*R2 ; 

0 

1 

-2*R 

3*R2 

-4*R*R2 

5*R2*R2; 

0 

0 

2 

-6*R 

12*R2 

-20*R*R2] 

« 

RMI-inv (RM) 
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EXACTSOL 


%  EXACTSOL. M  routine  for  PC-MATLAB 

% 

diary  exactsol . log 
( f ieldl , f ield2 , be ] =gmat lmat 
zero-zeros (6,6); 

* 

zvl-expmat ( tm ( 1 ) ) ; 
zv2 -expmat ( tm { 2 ) ) ; 
zv3=expmat  ( tm ( 3 ) ) ; 

« 

bcl=[zvl  zero  zero; 

eye ( 6 )  zero  zero ; 

zero  zv2  zero ; 

zero  eye ( 6 )  zero ; 

zero  zero  zv3 ; 

zero  zero  eye  ( 6 )  ] 

% 

bcO-bc*bcl 
CondNum-rcond ( bcO ) 
bcinv-inv(bcO) 

% 

force- [0;0;.5;0;0;0;0;0;0;0;0;0;0;0;0;0;0;.5] 

« 

UVW=bcinv* force 
uvw-zeros (7,11); 

uvw(7, 1; 11)*( . 5, .46, . 5-tm(l) , tm(2 ) /2 , .04, 0. , -.04,-tm(2) /2, tm(3)-. 5, -.46,-. 5 ] ; 
* 

UVWset- [ UVW (1:6) ] ; 
tmBottom-. 5-tm( 1 ) ; 

« 

for  ii-(lsll) 
if  ii— 4 

UVWset— ( UVW (7:12) ] ; 
tmBottom— tm  ( 2 )  /2 ; 
elseif  ii— 9 

UVWset- [UVW( 13: 18) ] ; 
tmBottom— .  5 ; 

end 

uvw ( 1 : 6 , ii ) -expmat ( uvw ( 7 , ii ) -tmBottom ) *UVWset ; 

end 

uvw 

« 

diary  exact sol. prn 
uvw' 

% 

stress-zeros (7, 11) ; 

stress (7,1: 11 )-[ .5, .46, .5-tm(l) , tm(2)/2, .04,0. , 

- . 04 , -tm (2)/2,tm(3)-.5,-.46,-.5J; 
strain-zeros (7,11); 

strain(7,l:ll)-[ .5, .46, .5-tm(l) ,tm(2)/2, .04,0. , 

-.04, -tm(2 )/2,tm(3)-.5,-.46,-.5); 

« 

StrainMat-(0  0  0  -p  00; 

000  0  -q  0; 

001  0  00; 

0  0  0  q  p  0; 

1  0  0  0  Op; 

010  0  0  q); 

« 

for  ii«(l:ll) 

strain ( 1 : 6 , ii ) -StrainMat  *uvw ( 1 : 6 , ii ) ; 
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stress ( 1 : 6 , ii ) =KM*  Str ainMat *uvw ( 1 s  6 , ii ) ; 

end 

diary  exstress. pm 

strain' 

stress' 

« 

diary  off 
if  iso“0 

uvwf -zeros (3,2); 

stress^zeros (7,2); 

stress (7,l:2)=[tf(l)  tf(2)J; 

strain-zeros (6,2); 

strain (6,l:2)=[tf(l)  tf(2)J; 

« 

StrainMat= [ -p  00; 

0  -q  0; 
q  p  0; 

0  0  p; 

0  0  q); 

« 

for  ii=(l:2) 

uvwf (l,ii) =UVW(ii* 6-2) ; 
uvwf (2,ii)=UVW(ii*6-l) ; 
uvwf (3,ii)=UVW(ii*6) ; 

strain ( 1 : 5 , ii ) =StrainMat *uvwf ( 1 : 3 , ii ) ; 

end 

stress (1:6,1) *SMl*StrainMat*uvwf (1:3,1); 

stress (1:6,2) =SM2  *  St rainMat  *uvwf (1:3,2); 

diary  ex-fstr.prn 

strain' 

stress' 

% 

diary  off 

end 


STXPL1ES 

%  STXPLIES  -  The  name  stands  for  "STRATA  THEORY ,  CROSS-PLIES".  This  routine 
%  creates  the  matrix  of  coefficients  that  arise  from  the  PDEs 

%  when  using  the  assumed  form  of  psi,  zeta,  and  w. 

% 

cl-matprops  ( 1 ) ; 
c 2 -mat props ( 2 ) ; 
c3»matprops ( 3 ) ; 
ps—p*p; 
qs-q*q; 

pq-p*q; 
bs-bet a*bet a ; 

stheory-zeros ( 2  *N+3 , 2  *N+3 ) ; 

stheory ( 1 , 1 ) * ( c3*ks ( 0 ) +ss ( 7 , 0 ) ) *ps+ ( c3  *ks ( 0 ) +ss (8,0)) *qs ; 

% 

for  ii«(l:N) 

stheory ( 1 , ii+2 ) «ii*beta*c3*p*ks ( ii-1 ) ; 
stheory ( 1 , ii+N+3 ) *ii*beta*c3*q*ks ( ii-1 ) ; 

end 

% 

for  jj«(0:N) 

stheory ( j  j  +2 , 1 ) * j  j *beta*c3*p*ks ( j  j-1 ) ; 
for  ii*(0:N) 

kst if f *ks ( ii+ j  j )  ; 

stheory ( j  j+2 , ii+2 ) =ii* j  j *bs*c3*ks ( ii+ j  j-2 ) 

+(ss(l,ii+j j )+cl*kstiff )*ps+(ss(4,ii+j j )+c3*kstiff )*qs; 
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Btheory ( j  j  +2 , ii+N+3 )  =  ( ss ( 3 , ii+ j  j ) +ss ( 4 , ii+ j j )  +  ( c2+c3 ) *kst if f ) *pq; 

end 

end 

% 

for  jj=(0:N) 

stheory ( j j+N+3, l)=j j*beta*c3*q*ks( j j-1) ; 
for  ii*(0:N) 

kstiff=ks(ii+jj) ; 

Btheory ( j  j  +N+3 , ii+2 ) = ( as ( 3 , ii+ j  j ) +sb ( 4 , ii+ j  j )  +  ( c2+c3 ) *kst if f ) *pq; 
stheory ( j  j  +N+3 , ii+N+3 ) =ii* j  j  *bs*c3*ks ( ii+ j  j -2 ) 

+(B8(2, ii+j j )+cl*kstiff )*qs+(ss(4, ii+j j )+c3*kstiff ) *ps; 

end 

end 

« 

diary  strtheor.log 
stheory 

f orce*zeros ( 2  *N+3 , 1 ) ; 
force(l)»l; 

CondNum*r cond ( Btheory ) 

% 

UVW“inv ( stheory ) *  force 
uvw=zeros (7,11); 

uvw(7,l:ll)  =  [ .5, .46, .5-tm(l) ,tm(2)/2, .04,0. ,-.04,-tm(2)/2,tm(3)-.5,-.46,-.5]  ; 

% 

alphafactor=( 1,1, 1,0, 0,0,0, 0,-1, -1,-1) ; 
for  ii-(lsll) 
uvw(l,ii)=0; 
uvw(2, ii)=0; 
uvw(4, ii)=0; 
uvw(5,ii)=0; 

apluabz=alpha*alphaf actor ( ii ) +bet a*uvw ( 7 , ii ) ; 
for  jj*(0:N) 

uvw(l,ii)=uvw(l, ii)+uvw(j j+2)*5j*beta*aplusbzA ( j j-1) ; 
uvw(2,  ii)=uvw(2,  ii)+UVW(j  j+N+3)  *  j  j*beta*aplusbz‘*  ( j  j-1) ; 
uvw(4,  ii)»uvw(4,  ii)+UVW(  j  j+2)  ♦apluBbz''  j  j ; 
uvw  ( 5 ,  ii )  =uvw  ( 5 ,  ii )  +UVW  ( j  j  +N+3 )  ♦apluBbz''  j  j  ; 

end 

uvw(3, ii)*0; 
uvw ( 6 , ii ) *UVW ( 1 ) ; 

end 

uvw 

diary  strtheor.prn 
uvw' 

« 

atresB=zeroB(7, 11) ; 

stress (7, 1:11)«[ .5, ,46,0.50-tm(l) ,tm(2)/2, .04,0. , 

-.04,-tm(2)/2,tm(3)-0.50,-.46,-.5J ; 
strain^zeros (7,11); 

strain(7,l:ll)*[ .5, ,46,0.50-tm(l) ,tm(2)/2, .04,0. , 

-.04,-tm(2) /2 , tm{3) -0. 50, -. 46,-. 5) ; 

% 

StrainMat*[0  0  0  -p  00; 

000  0  -q  0; 

0  0  1  0  0  0; 

0  0  0  q  p  0; 

1  0  0  0  Op; 

0  1  0  0  0  q); 

% 

for  ii-(lsll) 

strain ( 1 : 6 , ii ) «StrainMat *uvw ( 1 : 6 , i i ) ; 
stress ( 1 : 6, ii)“KM*StrainMat*uvw( 1:6, ii) ; 

end 

diary  ststress.prn 
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strain' 

stress' 

« 

diary  off 
if  iso==0 

uvwf “zeros (5,2); 

stress“zeros (7,2); 

stress (7,l:2)=(tf(l)  tf(2)]; 

strain=  zeros (6,2); 

strain (6,l:2)“[tf(l)  tf(2)); 

% 

StrainMat=[0  0  -p  00; 

00  0  -q  0; 

0  0  q  p  0; 

1C  0  Op; 

01  0  0  q]; 

« 

for  ii*(l:2) 

uvwf (l,ii)=0; 
uvwf  (2,ii)“0; 
uvwf (3,ii)=0; 
uvwf (4, ii j“0; 
zf“zdist(ii) ; 
for  jj“(0:N) 

uvwf  (1, ii)“uvwf (1, ii)+UVW( j j+2)* j j*zf A ( j j-1) ; 
uvwf (2,ii)=uvwf <2,ii)+UVW( jj+N+3)*jj*zfA( jj-1); 
uvwf ( 3 , ii ) =uvwf ( 3 , ii ) +UVW ( j  j  +2 ) *  zf  A  j  j ; 
uvwf (4, ii)*uvwf (4, ii)+UVW( j j+N+3) *zf A j j ; 

end 

uvwf (5,ii)=UVW(l); 

strain ( 1 : 5 , ii ) “StrainMat *uvwf ( 1 s  5 , ii ) ; 

end 

stress (1:6,1) =SMl*StrainMat *uvwf (1:5,1); 

stress (1:6,2) =SM2  *StrainMat *uvwf (1:5,2); 

diary  st-fstr.prn 

strain' 

stress ' 

t 

diary  off 

end 


EXP  MAT 

function  em=expmat ( z ) 

%  EXPMAT  -  Function  to  generate  the  matrix  exponential  for  the  boundary 
%  condition  matrices  of  the  "EXACT  STRATIFIED  PLATE  SOLUTION"  for 

%  the  analysis  of  composite  panels. 

ezp*exp(R*z) ; 
ezm»exp(-R*z) ; 

eVector*[ezp;  z*ezp;  z*z*ezp;  ezm;  z*ezm;  z*z*ezm]; 
aV«RMI *eVector 

em»aV ( 1 ) *eye ( 6 ) +aV( 2 ) *FM+aV ( 3 ) *FM*FM+aV ( 4 ) *FM*FM*FM+aV ( 5 ) *FM*FM*FM*FM 
+ aV ( 6 ) *  FM*  FM*  FH*  FM*  FM ; 


GMATLMAT 

f unct ion  [ f ml , f m2 , be ] »gmat lmat 

%  GMATLMAT  Function  to  generate  the  material  property  matrices  for  the 
%  "EXACT"  solution  of  the  cross-ply  composite  panel  problem. 

% 
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%  Syntax:  [ FMl , FM2 , BC ] =GMATLMAT 

« 


P 

* 

nm/a 

% 

q 

mtr/b 

« 

mp 

* 

Vector  of  material  properties  (1x3) 

« 

A 

fp 

SC 

Vector  of  fiber  strata  properties  (1x9) 

« 

BC 

s 

Boundary  condition  matrix 

« 

FMl 

* 

Field  matrix  #1 

% 

FM2 

s 

Field  matrix  #2 

t 

« 

where  ( FMl ] 

u  =  Displacement  vector 
u  =  L  (FM2 ]  u  L  =  Eigenvalue  (aka  lambda) 

% 

%  This  function  is  used  in  conjunction  with  other  MATLAB  routines  to 
%  verify  the  "STRATA  THEORY "  of  composite  panel  analysis.  Written  by 
%  Alan  L.  Lesmerises. 


cl«mat props ( 1 ) 
c2=mat props ( 2 ) 
c3»raatprops ( 3 ) 
c4*c2+c3 
z«zeros(3,3) ; 

« 

ml*[c3  00; 

0  c3  0; 

0  0  cl ] ; 

% 

m2*c4*[  0  0  p; 

0  0  q; 

-p  -q  0]; 

« 

m3*-[ (p*p*cl+q*q*c3)  p*q*c4 

p*q*c4  (p*p*c3+q*q*cl) 

0  0 

« 

%  Above  matrices  appear  in  PDE  as: 


0 

0  ; 
(p*p+q*q) *c3  j ; 


[ml]u"+[m2]u'+[m3]u=0 

or 


[ml]  [0] 

Lu 

[0] 

[ml] 

Lu 

[0]  -[m3] 

u 

=  L 

l  [ml] 

[m2]  j 

u 

. 

fml«[ml  z  ; 
z  -m3 ] ; 

% 

fm2» [  z  ml; 
ml  m2 ] ; 

% 

FH*inv( fm2 ) *fml 

% 

m2*[  0  0  p*c3; 

0  0  q*c3; 

-p*c2  -q*c2  0  ] ; 

« 

sl«f ibprops{ 1, 1) ; 
s2*f ibprops (1,2); 
s3»f ibprops (1,3); 
s4«  f ibprops (1,4); 
s7«f ibprops (1,7) ; 
s8*f ibprops (1,8); 

m2fl*m2+tf (1)*( (p*p*sl+q*q*s4)  p*q*(s3+s4)  0 
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% 


p*q* (b3+b4) 
0 


(p*p*s4+q*q*s2 ) 
0 


0 

(p*p*s7+q*q*s8) ) 


f 

9 


sl*f ibpropa (2,1); 
s2»f ibprops ( 2 , 2 ) ; 
s3~fibprops(2,3) ; 
s4=*f  ibprops  (2,4); 
s7«fibprops(2,7) ; 

■8*f ibprops (2, 8) ; 
m2f2*m2+tf (2)*[ (p*p*sl+q*q*s4) 
p*q* (s3+s4) 

0 


% 


p*q* ( s3+s4 ) 
(p*p*s4+q*q*s2 ) 
0 


0 

0 

(p*p*s7+q*q*s8) ) 


I3*eye(3) 
be* [ml  m2 
z  z 
z  z 
z  z 
z  z 
z  z 

% 


z  z 
-ml  -m2 
z  13 
z  z 
z  z 
z  z 


z  Z  Z  Z  z  z 

ml  m2fl  z  z  z  z 

z  -13  z  z  z  z 

z  z  -ml  -m2  ml  m2f2 
z  z  z  13  z  -13 

Z  Z  Z  Z  Z  z 


z  z; 
z  z; 
z  z; 
z  z; 
z  z; 
-ml  -m2) 


%  End  Of  function  GMATLMAT.M 

% 


TREND 

•etparam 

extrend*zeros (4,50); 
for  jndex*l:50 
A* jndex 
B«A; 

extrend ( 1 , jndex ) =A; 
extrend ( 2 , j ndex ) =B ; 
calparam 
exacteol 

extrend(3, jndex)=uvw(6,6)  ; 
atxplies 

extrend(4, jndex)*uvw(6, 6) ; 

end 

diary  trend. prn 
extrend ' 
diary  off 
quit 


KS 

function  k*ke( order) 

%  KS  -  Evaluates  the  integral  in  z  that  is  associated  with  the  k's  that 
%  appear  in  the  P.D.E.'s  for  the  "STRATA  THEORY"  for  the  analysis  of 

%  composite  panels. 

k-0; 

if  order  <  0 

elseif  (order/2  **  fix(order/2) \ 
k*l/ ( order+1 ) /2 "order /beta ; 

end 
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55 


function  8stiff=ss(sterm, order) 

t 

%  SS  -  Sum  the  terms  of  £  s(f,j)  *  z^n 

« 

%  Syntax  SS(j, order) 

% 

%  where: 


« 

j) 

sc 

membrane  elasticity  terns 

« 

f 

SC 

fiber  number 

% 

j 

= 

the  elasticity  tern  of  interest 

« 

a 

n 

s 

the  "order"  of  the  z  tern 

*0 

« 

stem 

= 

membrane  properties  vector 

« 

order 

= 

the  order  of  the  z  tern  (n  above) 

[n,numfibers]-size(zdist) ; 
sstif f~0; 

for  index=* ( 1 : numf ibers ) 

satif  f=sstif  f+f ibprops ( index,  Bterm)  *zdist  ( index)  ',order*tf  (index) ; 

end 

« 

%  End  of  SS.M  function 

« 
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Appendix  B:  MATLAB  Routines  for  Composite  Panel  in  Cylindrical  Bending 


SOLVE 

A«10 

aetparam 

calparam 

exactsol 

atxpliea 

quit 


SETPARAM 

%  SETPARAM. M  routine  for  PC-MATLAB 

« 

format  short  e 

global  n  A  p  mat props  fibprops  zdiat  N  alpha  beta  tf  tm  FM  RMI 
diary  setparam.log 

% 

iso=0 

n=l 

% 

Em=0.3e6 
NUm= . 40 

Gm-Em/( 1+NUm) /2 

« 

kl=Em* ( 1-NUm) / ( l-2*NUm) / { 1+NUm) 
k2“Em*NUm/ ( l-2*NUm) / ( 1+NUm) 
k3-Gm 

% 

KM«[kl  k2  0; 
k2  kl  0; 

0  0  k3] ; 

% 

matprops» [ kl  k2  k3] 

« 

if  iao»«l 
tf«[0  0] 
tm*[ .25  .5  .25] 
alpha*0 
beta-1 
« 

fibprops=[0  00000000; 

000000000]; 

elae 

aqrt2»aqrt ( 2 ) 

volf-aqrt2/2 

volm-l-volf 

« 

tf*(sqrt2/4  aqrt2/4] 

tm«[ .25-aqrt2/8  .5-sqrt2/4  .25-aqrt2/8] 

alpha—  ( l+aqrt2 )  /2 

beta-2+aqrt2 

« 

Ef-20e6 
HUf- . 35 

Gf«Ef/(l-*-NUf  )/2 

% 

al-volf *Ef * ( 1-NUf ) / { l-2*NUf ) / { 1+NUf ) +volm*kl 
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s2*kl/volm 

s3*volf *Ef *NUf / ( l-2*NUf ) / ( 1+NUf ) +volm*k2 
s4*Gni/volm 
■5*0 
■6*0 

■7*volf *Gf+volm*Gm 

■8*Gm/volm 

■9*0 

■10*k2 

% 

fibprops=[sl  s2  s3  s4  s5  s6  s7  b8  b9; 

b2  si  s3  b4  b5  s6  b8  s7  s9 ] 

% 


« 

SMI* [Si 

s3 

0 

0 

0; 

« 

S3 

s2 

0 

0 

0; 

% 

s3 

slO 

0 

0 

0; 

« 

0 

0 

84/2 

0 

0; 

« 

0 

0 

0 

s7/2 

0; 

« 

0 

0 

0 

0 

sS/2] 

« 

« 

SM2*(s2 

s3 

0 

0 

0; 

% 

s3 

si 

0 

0 

0; 

« 

slO 

S3 

0 

0 

0; 

« 

0 

0 

84/2 

0 

0; 

% 

0 

0 

0  - 

s8/2 

0; 

% 

0 

0 

0 

0  - 

s7/2] 

« 

SMI* [si  0; 

0  s7/2 ] 

« 

SM2* [ b2  0; 

0  b8/2  ] 

end 

% 

zdist*[ .25,  -.25] 
N-5 


CAL  P ARAM 

*  C ALP ARAM. M  routine  for  PC-MATLAB 

« 

p*n*pi/A 

ps-p*p 


RM*[  1 

P 

ps 

p*ps; 

0 

1 

2*p 

3*ps; 

1 

-P 

ps 

-p*ps; 

0 

1 

-2*p 

3*ps] 

RMI*inv(RM) 


EXACTS OL 

%  EXACTSOL . M  routine  for  PC-MATLAB 

« 

diary  exactsol.log 
[ f ieldl , f ield2 , be ] *gmat lmat 
zero*zeros (4,4); 

« 

zvl*expraat ( tm ( 1 ) ) ; 
zv2*expmat ( tro ( 2 ) ) ; 
zv3*expmat ( tm ( 3 ) ) ; 

« 
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bcl*[zvl  zero  zero; 

eye ( 4 )  zero  zero ; 

zero  zv2  zero; 

zero  eye ( 4 )  zero ; 

zero  zero  zv3 ; 

zero  zero  eye ( 4 ) ] 

% 

bc0=bc*bcl 
CondNum-r cond { bcO ) 
bcinv-inv ( bcO ) 

% 

force**  [0;  .5;0;0;0;0;0;0;0;0;0;.5] 

% 

UW=bcinv*force 
uw- zeros (5, 11) ; 

uw(5,l:ll)=( .5, .46, .5-tm(l) ,tm(2)/2, .04,0. ,-.04,-tm(2)/2,tm(3)-.5,-.46, 

% 

UWset- [ UW ( 1 : 4 ) ] ; 
tmBottom= - 5-tm ( 1 ) ; 

% 

for  index-(l:  11) 
if  index— 4 

UWset= [ UW (5:8) ] ; 
tmBott om=-tm ( 2 ) /2 ; 
eleeif  index==9 

UWset *=  [  UW  ( 9 : 12 )  ]  ; 
tmBottom=- . 5 ; 

end 

uw ( 1 : 4 , index ) =expmat ( uw ( 5 , index ) -tmBottom ) *UWset ; 

end 

uw 

« 

diary  exactsol.prn 
uw' 

diary  off 

% 

stress=zeros(4, 11) ; 

stress (4, 1: 11) =[ .5, .46, .5-tm(l) ,tm(2)/2, .04,0. , 

-.04,-tm(2)/2,tm(3)-.5,-.46,-.5] ; 

8train=*zero8(4, 11) ; 

strain(4,l:ll)*( .5, .46, .5-tm(l) ,tm(2)/2, .04,0. , 

-.04,-tm(2)/2,tm(3)-.5,-.46,-.5] ; 

StrainMat=[0  0  -p  0; 

0  1  0  0; 

10  Op]; 

% 

for  ii-(lsll) 

strain(l:3, ii)*StrainMat*uw(l:4, ii) ; 
stress ( 1 : 3 , ii ) *KM*StrainMat *uw ( 1 : 4 , ii ) ; 

end 

diary  exstress. pm 

strain' 

stress' 

diary  off 

« 

if  iso— 0 

uwf -zeros ( 2 , 2 ) ; 
stress-zeros (3,2); 

•tress (3, l:2)-[tf (1)  tf(2)]; 
strain-zeros (3,2); 

•train(3,l:2)«[tf (1)  tf(2)]; 

* 
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StrainMat* [ -p  0; 

0  p]; 

« 

for  iis(l:2) 

uwf (l,ii)=UW(ii*4-l) ; 
uwf { 2 , ii ) =UW ( ii*4 ) ; 

strain ( 1 : 2 , ii ) =StrainMat*uwf ( 1 : 2 , ii ) ; 

end 

stress (1:2 , l)=SMl*StrainMat*uwf (1:2, 1) ; 

stress (1:2,2) =SM2*StrainMat*uwf (1:2,2); 

diary  ex-fstr.prn 

strain' 

stress' 

« 

diary  off 

end 


STXPL1ES 

*  STXPLIES  -  The  name  stands  for  "STRATA  THEORY,  CROSS-PLIES".  This  routine 
%  creates  the  matrix  of  coefficients  that  arise  from  the  PDEs 

%  when  using  the  assumed  form  of  psi  and  w. 

% 

cl»mat props ( 1 ) ; 

c2»matprops ( 2 ) ; 

c3=matprops ( 3 )  ; 

ps-p*p; 

bs«beta*bet a ; 

stheory*zeroB(N+2,N+2) ; 

stheory  ( 1 , 1  *  ■«  ( c3*ks  ( 0 )  +ss  (7,0))  *ps ; 

« 

for  ii-(l:N, 

stheory  (1  *-2  )«ii*beta*c3*p*ks(ii-l) ; 

end 

% 

for  jj*(0:N) 

stheory ( j  j  +2 , 1 ) * j  j *beta*c3*p*ks ( j  j— 1 ) ; 
for  ii-(0:N) 

stheory ( j j+2, ii+2)=ii* j j*bs*c3*ks(ii+j j-2)+(ss(l, ii+j j )+cl*ks(ii+j j ) ) *ps; 

end 

end 

« 

diary  strtheor.log 
stheory 

f orce* zeros (N+2 , 1 ) ; 
force(l>«l; 

CondNum*rcond ( stheory ) 

UW»inv ( stheory ) *  force 
uw*zeros(5, 11) ; 

uw(5,l:ll)*[ .5, .46, .5-tra(l) ,tm(2)/2, .04,0. ,-.04,-tm(2)/2,tm(3)-.5,-.46,-.5] ; 
alphaf actor* ( 1, 1, 1,0, 0,0,0, 0,-1, -1,-1] ; 
for  ii-(ltll) 
uw(l,ii)*0; 
uw(3, ii)*0; 

aplusbz»alpha*alphaf actor (ii)+beta*uw(5,  ii) ; 
for  jj*(0:N) 

uw( 1, ii)*uw(l, ii)+UW{ j j+2 ) * j j*beta*aplusbzA ( j j-1) ; 
uw( 3 , ii ) *uw ( 3 , ii ) +UW ( j  j  +2 ) *aplusbz A  j  j ; 

end 

uw(2, ii)*0; 
uw(4,ii)«UW(l) ; 

end 
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uw 

diary  strtheor.prn 
uw' 

diary  off 

% 

stress* zeros ( 4 , 11) ; 

stress (4, 1: 11)=( . 5, .46,0. 50-tm( 1) , tm(2 ) /2 , .04, 0. , 

-.04,-tm(2)/2,tm(3)-0.50,-.46,-.5] ; 
atrain=zeros (4, 11 ) ; 

strain ( 4, 1: 11)  =  [ .5, .46,0.50-tm(l) , tm(2)/2, .04,0. , 

-.04,-tm(2)/2,tm(3)-0.50,-.46,-.5] ; 

« 

StrainMat=[0  0  -p  0; 

0  1  0  0; 

10  0  p); 

% 

for  ii-(lsll) 

strain ( 1 : 3 , ii ) =StrainMat  *uw ( 1 : 4 , ii ) ; 
stress (1:3, ii)=KM*StrainMat*uw(l:4, ii) ; 

end 

diary  ststress.prn 

strain' 

stress' 

% 

diary  off 
if  iso*=0 

uwf =zeros (3,2); 
stress=zeros (4,2); 
stress(4,l:2)=[tf (1)  tf(2)]; 
strain=zeros (4,2); 
strain(4, 1:2 )=(tf ( 1)  tf (2 ) ] ; 

« 

StrainMat=[0  -p  0; 

1  Op); 

« 

for  ii=(l:2) 
uwf (1, ii)=0; 
uwf (2 , ii)=0; 
zf =zdist ( ii ) ; 
for  jj=(0:N) 

uwf (1, ii)=uwf (1, ii)+uw( j j+2)* j j*zf A ( j j-1) ; 
uwf (2,ii)=uwf (2,ii)+UW{ j j+2)*zf A j j; 

end 

uwf ( 3 , ii ) — UW ( 1 ) ; 

strain ( 1:2 , ii)-StrainMat*uwf ( 1: 3 , ii) ; 

end 

stress (1:2,1) =SMl*StrainMat*uwf (1:3,1); 
stress (1:2,2) *SM2*StrainMat*uwf (1:3,2); 
diary  st-fstr.prn 
strain' 
stress ' 

% 

diary  off 

end 


EXPMAT 

function  em=expmat ( z ) 

%  EXPMAT  -  Function  to  generate  the  matrix  exponential  for  the  boundary 
%  condition  matrices  of  the  "EXACT  STRATIFIED  PLATE  SOLUTION"  for 

%  the  analysis  of  composite  panels. 
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ezp*exp(p*z) ; 
ezm»exp(-p*z) ; 

eVector=[ezp;  z*ezp;  ezm;  z*ezm] ; 
aV=RHI *eVector 

em=aV ( 1 ) *eye ( 4 ) +aV ( 2 ) *FM+aV( 3 ) *FM*FM+aV( 4 ) *FM*FM*FM; 


GMATLMAT 


function  [fml,fm2,bc]=gmatlmat 

%  GMATLMAT  Function  to  generate  the  material  property  matrices  for  the 
%  "EXACT"  solution  of  the  cross-ply  composite  panel  problem. 


%  Syntax:  [ FM1 , FM2 , BC ] =GMATLMAT 


« 

« 

% 

% 

« 

« 

« 

% 

% 

« 

% 


p 

= 

nm/a 

mp 

= 

Vector  of  material  properties  (1x3) 

fp 

= 

Vector  of  fiber  strata  properties  (1x9) 

BC 

= 

Boundary  condition  matrix 

FM1 

= 

Field  matrix  #1 

FM2 

= 

Field  matrix  #2 

u  =  Displacement  vector 

where  [ FM1 ]  u  =  L  [FM2 ]  u  L  =  Eigenvalue  (aka  lambda) 


%  This  function  is  used  in  conjunction  with  other  MATLAB  routines  to 
%  verify  the  "STRATA  THEORY"  of  composite  panel  analysis.  Written  by 
%  Alan  L.  Lesmerises. 


cl=matprops ( 1 ) ; 
c2=matprops (2 ) ; 
c3=matprops ( 3 ) ; 
c4*=c2+c3 ; 
z^zeros (2,2) ; 

% 

ml=[c3  0; 

0  cl]; 

% 

m2=c4*p*[0  1; 

-1  0]; 

% 

m3=(-p*p*cl  0  ; 

0  -p*p*c3); 

% 

%  Above  matrices  appear  in  PDE  as: 

% 

%  [ml)u"+(m2]u'+[m3]u=0 

*  or 


[ml]  [0] 

Lu 

[0] 

[ml] 

Lu 

[0]  -[m3] 

u 

*  L 

(ml] 

[m2] 

u 

% 

fml=[ml  z  ; 
z  -m3 ) ; 

% 

fm2*[  z  ml; 
ml  m2 ] ; 

% 

FM*inv(fm2)*fml 

% 

ml*[c3  0; 

0  cl]; 
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« 

m2=[  0  p*c3; 

-p*c2  0  ] ; 

* 

Bl=f ibprops (l,l)*tf (1) ; 
s2*fibprops(l,2)*tf (1) ; 
s3=f ibprops ( 1 , 3 ) *t f ( 1 ) ; 
e4=f ibprops ( 1 , 4 ) *t f ( 1 ) ; 
s7=fibprops(l,7)*tf (1) ; 
a8=f ibprops (1,8) *tf (1) ; 
m2f l“m2+(p*p*sl  0  ; 


Bl*f ibprops ( 2 , 1 ) *t f ( 2 ) ; 
a2=fibprops(2,2)*tf (2) ; 
b 3= f ibprops ( 2 , 3 ) *tf ( 2 ) ; 
s4*fibprops(2,4)*tf (2) ; 
a7*fibprops(2,7)*tf (2) ; 
s8*fibprops(2,8)*tf (2) ; 
m2f2*m2+[p*p*8l  0  ; 

0  p*p*s7] ; 

% 


I2=eye(2) ; 


be- [ml  m2 

z 

z 

z 

z 

z 

z 

z 

z 

z 

z; 

z  z 

-ml 

-m2 

ml 

m2f  1 

z 

z 

z 

z 

z 

z; 

z  z 

z 

12 

z 

-12 

z 

z 

z 

z 

z 

z; 

z  z 

z 

z 

z 

z 

-ml 

-m2 

ml 

m2f  2 

z 

z; 

z  z 

z 

z 

z 

z 

z 

12 

z 

-12 

z 

z; 

z  z 

« 

z 

z 

z 

z 

z 

z 

z 

z 

-ml 

-m2  ] 

%  End  of  function  GMATLMAT.M 
% 


TREND 

setparam 

extrend*zeros (3,50); 
for  jndex«l:50 
A* jndex 

extrend ( 1 , j ndex ) =A ; 

calparam 

exactsol 

extrend(2, jndex)=uw(4,5) ; 
Btxpliee 

extrend  ( 3 ,  j  ndex )  *=uw  (4,5); 

end 

diary  trend.prn 
extrend ' 
diary  off 
quit 


KS 

function  k«ks( order) 

%  KS  -  Evaluates  the  integral  in  z  that  is  associated  with  the  k's  that 
%  appear  in  the  P.D.E.'s  for  the  "STRATA  THEORY"  for  the  analysis  of 

%  composite  panels. 


k«0; 

if  order  <  0 

elseif  (order/2  fix ( order /2) ) 
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k=l/ (order+1) /2*order/beta; 

end 


SS 

function  satiff=ss(stenn/ order) 

% 

%  SS  -  Sum  the  terms  of  £  s(f,j)  *  zAn 

% 

*  Syntax  SS(j, order) 

% 

%  where: 

%  s(f,j)  =  membrane  elasticity  terms 

%  f  =  fiber  number 

%  j  =  the  elasticity  term  of  interest 

%  n  =  the  "order"  of  the  z  term 

« 

%  stem  =  membrane  properties  vector 

%  order  =  the  order  of  the  z  tern  (n  above) 

[n,numfibers]=size(zdist) ; 
sstiff=0; 
s sum-0; 

for  index* (l:numfibers) 

sstif  f =sst if  f +f  ibprops  ( index ,  stem)  *zdist  ( index )  *order*tf  ( index ) ; 

end 

« 

%  End  of  SS.M  function 

« 


98 


Bibliography 


1.  Chen,  Chi-Tsong,  Introduction  to  Linear  System  Theory.  New  York:  Holt, 
Rinehart  and  Winston,  Inc.,  1970. 

2.  Jones,  Robert  M.,  Mechanics  of  Composite  Materials.  New  York:  Hemisphere 
Publishing  Corp.,  1975. 

3.  Kapania,  Rakesh  K.  and  Stefano  Raciti.  "Recent  Advances  in  Analysis  of 
Laminated  Beams  and  Plates,  Part  II:  Vibrations  and  Wave  Propagation,"  AIAA 
Journal,  27:  935-946  (July  1989). 

4.  Kulkami,  S.V.  and  N.J.  Pagano.  "Dynamic  Characteristics  of  Composite 
Laminates,"  Journal  of  Sound  and  Vibration,  23:  127-143  (July  1972). 

5.  Meirovitch,  Leonard,  Analytical  Methods  in  Vibrations.  New  York:  Macmillan 
Publishing  Co.,  1967. 

6.  Meirovitch,  Leonard,  Elements  of  Vibration  Analysis  (Second  Edition).  New 
York:  McGraw-Hill  Book  Company,  1986. 

7.  Mindlin,  R.D.  "Influence  of  Rotatory  Inertia  and  Shear  on  Flexural  Motions  of 
Isotropic,  Elastic  Plates,"  Journal  of  Applied  Mechanics,  18:  31-38  (March  1951). 

8.  Pagano,  N.J.  "Exact  Solutions  for  Composite  Laminates  in  Cylindrical  Bending," 
Journal  of  Composites,  3:  398-411  (July  1969). 

9.  Pagano,  N.J.  "Exact  Solutions  for  Rectangular  Bidirectional  Composites  and 
Sandwich  Plates,"  Journal  of  Composites,  4:  20-34  (January  1970). 

10.  Reissner,  Eric.  "The  Effect  of  Transverse  Shear  Deformation  on  the  Bending  of 
Elastic  Plates,"  Journal  of  Applied  Mechanics,  12:  A69-A77  (June  1945). 

11.  Saada,  Adel  S.,  Elasticity  -  Theory  and  Applications.  Malabar,  Florida:  Robert 
E.  Krieger  Publishing  Co.,  1974. 

12.  Whitney,  J.M.  and  N.J.  Pagano.  "Shear  Deformation  in  Heterogeneous 
Anisotropic  Plates,"  Journal  of  Applied  Mechanics,  37:  1031-1036  (December 
1970). 


99 


Vita 

Captain  Alan  L.  Les menses  was  bom  on  January  7,  1959  in  Manchester,  NH. 
He  graduated  from  Trinity  High  School,  Manchester,  NH,  in  1976.  He  then  attended 
Daniel  Webster  College  (1976-77),  majoring  in  Aeronautical  Engineering.  He  enlisted 
in  the  U.S.  Air  Force  in  December  1977,  and  was  assigned  to  Plattsburg  AFB,  NY  as 
a  Jet  Engine  Technician.  He  continued  his  education  at  Plattsburg  State  University 
College,  SUNY  (1977-79).  In  1983,  he  was  assigned  to  the  University  of  Oklahoma  after 
being  accepted  into  the  Airman  Education  and  Commissioning  Program.  He  graduated 
with  honors  in  May  1986  with  a  Bachelor  of  Science  degree  in  Aerospace  Engineering. 
After  attending  Officer  Training  School,  he  was  assigned  to  the  Fuels  and  Lubrication 
Division,  Aero  Propulsion  Laboratory,  Air  Force  Wright  Aeronautical  Laboratories, 
Wright-Patterson  AFB,  OH.  During  this  assignment,  he  was  program  manager  for 
several  research  and  development  contracts,  and  developed  computer  software.  He  also 
analyzed  test  data,  and  designed  test  equipment  and  facilities  for  fundamental  combustion 
research.  While  there,  he  co-authored  six  technical  papers  based  on  this  research. 

Permanent  Address: 

Alan  L.  Lesmerises 

62  Harrington  Ave. 

Concord,  MA  01742 


100 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-01 B8 


Publk  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  e**tmg  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services.  Directorate  for  information  Operations  a  no  heports,  1215  Jefferson 
Davis  Highway.  Suite  1 204,  Arlington.  VA  22202-4302.  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0 188).  Washington,  DC  20S03 


1.  AGENCY  USE  ONLY  (Leave  blank) 


2.  REPORT  DATE 

March  1992 


3.  REPORT  TYPE  AND  OATES  COVERED 

Master's  Thesis 


4.  TITLE  AND  SUBTITLE 

DEVELOPMENT  OF  A  STRATIFIED  PLATE  MODEL  FOR 
COMPOSITE  PANELS 


6.  AUTHOR(S) 

Alan  L.  Lesmerises,  Captain,  USAF 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADORESS(ES) 

Air  Force  Institute  of  Technology,  WPAFB,  OH,-  45433-6583 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

AFIT/GA/ENY/92M-02 


9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  AODRESS(ES) 

Mr.  Jerome  Pearson,  WL/FIBGD,  WPAFB  OH  45433-6563 
Dr.  T.  Nicholas,  WL/MLUSf,  WPAFB  OH  45433-6563 


10.  SPONSORING  /MONITORING 
AGENCY  REPORT  NUMBER 


12a.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 


12b.  DISTRIBUTION  CODE 


13.  ABSTRACT  (Maximum  200  words) 

A  new  motifeL  for  composites  represents  a  composite  panel  as  an  alternating  stack 
of  isotropic  matrix  and  orthotropic  fiber-dominated  layers.  An  energy  approach 
results  in  a  set  of  partial  differential  equations  and  boundary  conditions  where 
the  mechanical  properties  of  the  matrix  and  fibeEs  appear  separately.  Two 
solution  toethods  ace  developed.  The  first,  called  the  "Strata  Theory"  is  developed 
for  general  applications,  while  the  second,  called  the  "Exact  Stratified  Plate", 
is  used  as  a  benctaaafck  for  comparison  to  Strata  Theory.  Both  solutions  assume  zero 
gradients  of  displacement  through  the  thickness  of  a  fiber  layer.  They  are 
compared  for  an  isotropic  material  and  a  two-ply  composite  with  cross-plies,  using 
simply-supported  cases  of  a  square  plate  and  cylindrical  bending. 


14.  SUBJECT  TERMS 


Stratified  Plate,  Composite,  Composite  Model,  Cross-Plies 


15.  NUMBER  OF  PAGES 

111 


16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 

Unclassified 


18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

Unclassified 


19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

Unclassified 


20.  LIMITATION  OF  ABSTRACT 


NSN  7540-01-280-5500 


Standard  Form  298  (Rev  2-89) 

^scribed  by  ansi  239- '8 

-^07 


